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A Method for Simplify 


ing the Calculations of 


the Natural Frequencies for a System 
Consisting of n Rigid Rotating Discs 
Mounted on an Elastic Shaft? 


JOHN L. BOGDANOFFt 
Wright Aeronautical Corporation 


SUMMARY 


A. Stodola! has shown how the positive and negative syn- 
chronous natural frequencies (N.F.’s) for a system consisting of 
nrigid discs mounted on a rotating elastic shaft can be calculated. 
(The positive and negative synchronous N.F.’s are defined as the 
N.F.’s that coincide with or are exactly opposite to, respectively, 
the rotational speed of the shaft.) The principal assumption 
used in his method is that the elastic characteristics of the system 
are symmetrical. This assumption, combined with a need for only 
the synchronous N.F.’s, enabled Stodola to reduce the problem 
from that of solving for all the roots of an algebraic polynomial 
(in one unknown) of 8n degree to that of obtaining all the roots of 
a similar polynomial of 2n degree. 

L. Foppl,? using the same procedure as Stodola, obtained the 
positive and negative synchronous N.F.’s for a rigid disc mounted 
on a rotating cantilever shaft. In addition, he obtained a fourth 
degree polynomial (with coefficients that are a function of the 
disc velocity) whose roots give the N.F.’s of the system. How- 
ever, this polynomial was not solved explicitly. It was used to 
obtain limiting expressions for the N.F.’s when the disc angular 
velocity was large. 

In the present paper, the writer discusses a method for finding 
all the N.F.’s for a system consisting of m rigid rotating discs 
mounted on an elastic shaft. With the assumption of elastic 
symmetry, the problem is reduced to finding all the roots of a 
polynomial of 2n degree whose coefficients are functions of the 
angular velocities of the discs. It is also possible to construct an 
electronic network and a stationary mechanical model which give 
all and some, respectively, of the N.F.’s of the original (symmet- 
ric) system. If the system does not possess elastic symmetry, 
the method reduces the problem to that of finding all the roots of 
asimilar polynomial of 4n degree. For this case (unsymmetrical) 
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it is not possible to construct a stationary mechanical model. 
However, the electronic network can still be made, although it 
is more complicated than for the symmetric case. 


INTRODUCTION 


r THE PAST, the principal source of vibration excita- 
tion in a system consisting of m rigid rotating discs 
mounted on an elastic shaft was unbalance in the discs. 
As a result, the positive and negative synchronous 
N.F.’s were of primary interest. With the advent of 
the aircraft propeller-gas turbine combination, addi- 
tional exciting forces, generated by the propeller, act 
on the turbine rotating system. These forces can be re- 
solved, by harmonic analysis, into a Fourier series whose 
terms have periods that are some rational multiple 
(positive or negative) of one of the disc angular veloci- 
ties. Therefore, in order to design a system of this 
type which is free from resonance vibration, it is neces- 
sary to know all the natural frequencies of the system as 
a function of a disc velocity. 

If the radii of gyration of each disc are small in com- 
parison with the length of the shaft, the gyroscopic 
moments may be neglected and each disc may be con- 
sidered as a point mass. In this case, many methods 
are available for determining the m natural frequencies. 
Essentially, these methods consist in finding all the 
roots of a polynomial of mth degree, if the system pos- 
sesses symmetrical elastic properties. However, if the 


gyroscopic moments of each disc cannot be neglected, 
the only published procedures known to the writer for 
finding the N.F.’s of the system involve finding all the 
roots of a polynomial of 8 or 4n degrees whose coef- 
ficients are functions of the angular velocities of the 


— 


— 
AS 
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discs. (The polynomial is of 8n degree if the system has 
unsymmetrical elastic properties and 4m degree if these 
properties are symmetrical.) It is easily seen that in 
order to find all the N.F.’s, as a function of a disc veloc- 
ity, for a system with n 2 3, a great deal of labor must 
be expended. 

This paper presents a method whereby all the N.F.’s 
may be obtained by solving for all the roots of a similar 
polynomial of 4n or 2n degree (the polynomial is of 4n 
degree for a system with unsymmetrical elastic proper- 
ties and 2m degree if elastic symmetry exists). The 
first part applies this method to the determination of 
all the N.F.’s for a system consisting of a single rigid 
rotating disc mounted on a symmetrical elastic canti- 
lever shaft as a function of the disc angular velocity. 
The second part of the paper applies this method to the 
determination of all the N.F.’s for rigid rotating disc 
mounted on an elastic shaft as a function of the angular 
velocity of one of the discs. An Appendix has been 
included which contains a rigorous analysis of the pro- 
cedures used in the first and second parts. 


NOTATION 
Symbol Unit Significance 
c.g. aida Center of gravity of disc 
1 in. Length 
M, Mi, Ms, .... Ib.in.—1sec.2, Mass 
A, C, As, C2, As, Cs, .... Ib.in.sec.2 Mass moment of inertia 
Jat» in. Coordinates of c.g. 
disc axis 
B Ib.in.? Shaft flexibility 
@, W2, Wh, rad.sec.~! Angular velocities 
T, T2) rad.sec.~! Angular velocities 
o, G1, 02, rad.sec.~! Angular velocities 
2.718281284 


-1 


in.lbs Kinetic energy 

in.lbs Potential energy 

To, To’, Ti in.Ibs.sec.2, Reduced kinetic energy 
0, Vo’ in.Ibs. Reduced potential energy 

m,k,n Real numbers 

in.’?]b.12 Normal coordinates 

j Specific index 

r,s Running indexes 


I. Rucip Disc on EvAstic CANTILEVER SHAFT 


Fig. 1 shows a rigid rotating disc, D, mounted on a 
vertical cantilever shaft. The disc is mounted rigidly 
on the shaft, and it is assumed that the disc c.g., O’, 
and axis coincide with the elastic axis of the shaft when 
the system is at rest (as shown in Fig. 1). It is assumed 
that the cross section of the shaft is uniform along its 
length, /, and that the principal area moments of inertia 
of the cross section are equal. w is the constant angular 
velocity of the disc. Fig. 2 shows an instantaneous 
position of the system. The X YZ axes, with origin O 
at the built-in end of the cantilever, are fixed in space; 
OZ is vertical. The X’Y’Z’ axes, with origin O’ at the 
c.g. of the disc, remain parallel to the X YZ axes during 
motion. gi and gs; are the coordinates of O’ pro- 
jected on the XOY plane; gq, represents gu + 1@21. 
qu and gs are the coordinates of the projection of P’ 
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Fic. 1. Rigid rotating disc mounted on an elastic shaft. 


on the X’O’Y’ plane; q is used to represent gw + 
ig. (It should be noted that the length of O’P’ is 
taken to be one unit.) qu, ga, and giz, go2 are assumed 
to be small. MM is the mass of the disc; Cis the mo- 
ment of inertia about its axis; A is the remaining princi- 
pal moment of inertia. The shaft is assumed massless 


and B denotes its flexural rigidity. 
The kinetic and potential energies of the system 


are 


2T = M(gu? + Ga?) + A(di2” + Gao?) — 
— + Cw? (1) 


and 
12B 6B 
(qu? + ga”) — 2 + + 
4B 
(q12” + 22”) (2) 


respectively. Lagrange’s equations for this case are 
d ( or) oV 


dt Ogn 
dt Og;2 Ogr2 

(a1) 

dt\Ow 


The equations of motion, obtained by substituting 
Eqs. (1) and (2) in Eqs. (a), have the following form: 


Fic. 2. 


Examit 
shows 
fourth 


and 
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(1) 


(2) 
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Fic. 2. Instantaneous position of rigid rotating disc mounted 
on a cantilever shaft. 


Mgu = qu + 

Mgu = — gu + 
Agu + = qu — (3) 
Agua — = gu — 


Examination of the first four equations of Eqs. (3) 
shows that the first and second and the third and 
fourth may be combined as 

12B 6B 


6B 4B (4) 
Ag — = — qe 


Mi = 


q@ and gq are vectors in the XOY and X’O’Y’ planes, 
respectively. 

The first of Eqs. (4) is the vector form of the linear 
momentum equation. It shows that the motion is 


essentially two-dimensional. The second is the vector 


form for the two components of the moment of momen- 
tum equations that are in the X’O’Y’ plane. The re- 
mainder of the moment of momentum equation is given 
by: 

dt 
This is the last of Eqs. (3). It shows that the compon- 
ent of the angular momentum vector in the O’P’ (or 


O’Z') direction is constant. It becomes, after carrying 
out the indicated differentiation, 


| aun) | = 0 (b) 


— = O (5) 


Eqs. (4) are linear differential equations with con- 
stant coefficients. Assume a solution of the form 


irt 
qa = 
qe (6) 
Substituting these Eqs. (4), 
12B B 
(7) 
6B 4B 
—Ar*qu + = qi — 
The characteristic equation is 
j12B _ 6B 
Mr? 
=0 (8) 


Cwr — Ar? 


It will be shown below that this equation has four real 
roots (two positive and two negative) for any given 
value of w. Therefore, the solutions of Eqs. (4) are 
periodic. 

The N.F.’s for this system are given by the roots of 
Eq. (8). In expanded form, this equation is 


r? a? fr? TW 


where a? = 3B/MI', 8? = 2B/Al, p = C/A. 
The corresponding natural médes can be obtained from 
either 


= */2[4 — (7?/a?)] (10a) 


\a? a? 


The three parameters, a, a/8, and p, describe the physi- 


(10b) 


. cal characteristics of the system. They have the follow- 


ing meaning: 

(1) q@ is the natural frequency of the system when 
the disc, D, is assumed to be a point mass. 

(2) a/8 furnishes the relation between the length of 
shaft and the radius of gyration of the disc correspond- 
ing to the polar moment of inertia C. 
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(3) p describes the disc form. 
a and a/8 are both positive and are unrestricted in 
value. », however, must satisfy the relation: 


Eq. (6) shows that at a natural frequency a plane 
through OZ and O’ contains O’P’. Further, the plane 
OZO’ rotates with constant angular velocity 7 about 
the line OZ. Using these results, it is possible to derive 
Eqs. (7) by considering the equilibrium of the inertia 
and elastic forces in, and moments perpendicular to, 
the plane OZO’. This procedure was followed by 
Stodola! and Foppl,? the latter arriving at Eq. (9) for 
the N.F.’s. 

Eq. (9) has four real roots, two positive and two nega- 
tive, for a fixed but arbitrary value of w. If. these 
N F.’s, 7, are plotted as the ordinate with w as the 
abscissa, four continuous curves are obtained. This 
process of N.F. determination for a given system having 
fixed values for p, a?/8? is equivalent to finding all the 
common values of the following two functions: 


w/a = constant 


where the constant assumes, in turn, all values in the 
interval and +”, Geometrically, this means 
that the curves in the (w/a) — (r/a) plane, defined by 
the first of Eqs. (11), are found by obtaining their inter- 
sections with the line w/a = constant. This line, by 
sweeping out the entire (w/a) — (r/a) plane, may be 
considered as a scanning device. It is obviously im- 
material as to just what curve or line is used as the 
scanning device. The only requirement is that it can 
be made to pass through every point in the (w/a) — 
(r/a) plane. Therefore, replace Eqs. (11) with 


w/a = m(r/a) 


where m is a real number which varies from — © to 


+o, If the second of Eqs. (12) is substituted in the 


first, we obtain 
an 
a? a? 


Eqs. (13) and (10a), with a?/8? = 1 and p = 2, were 
used to obtain Figs. 3a and 3b, respectively. Fig. 
3a gives the N.F.’s, and Fig. 3b gives the corresponding 
modes of motion. . 

Actually, Figs. 3a and 3b are only the upper half of 
the complete diagram for the N.F.’s and the modes of 
motion. However, all information can be obtained from 
these figures in the following manner. Let w, > 0 be 
the angular velocity of the shaft. Find the N.F.’s 
and modes corresponding to +, and —,; let these 


be denoted by 71, 72, 73, 7, and a, a2, as and ay, respec- 
tively. Then the N.F.’s and modes corresponding to 
this shaft speed are 71, 72, —73, —74 and a4, a2, a3, a. 
The case w; < 0 needs no comment, since the shaft 
speed may always be assumed positive. The lower 
and upper curves in Fig. 3a will hereafter be called the 
first and second mode of motion, respectively. A 
method will now be described for using Fig. 3 to de- 
termine the resonance conditions of the system. 

The most important exciting forces which can act 
on the system shown in Fig. 2 may be classified as 
follows: (1) periodic force acting in one plane; (2) 
equal periodic forces acting in two mutually perpen- 
dicular planes and out of phase by 7/2 rad.; (3) fre- 
quency of exciting force is independent of disc angular 
velocity and of type 1 or 2; (4) frequency of exciting 
force is related to disc angular velocity by a constant 
factor and of type 1 or 2. 

Analysis shows! that the distinction between (1) and 
(2) is not important for this case. Hence, an exciting 
force of type 2 will be assumed in what follows. 

If the frequency of the exciting force is of type 3, the 
disc angular velocities at which resonance will occur are 
found as follows: 


(a) In Fig. 3a, draw a line parallel to abscissa with 
ordinate equal to the exciting frequency divided by a. 

(b) Find the abscissa of the intersection of this 
line with the N.F. curves. 

(c) The abscissa of these intersections gives the 
dimensionless disc speeds, w/a, at which resonance 
will occur. 

(d) The mode corresponding to this N.F. and disc 
speed can be found directly from Fig. 3b. 

If the frequency of the exciting force is of type 4, 
the N.F.’s and corresponding disc speeds may be found 
as follows: 
(a) Assume that the frequency of the exciting force, 
a, is related to the angular velocity, w, of the disc by 
o=k,w. 
(b) In Fig. 3a, draw line ke(w/a) = 7/a. 
(c) Find the abscissa and ordinate of the inter- 
sections of this line with the N.F. curves. 
(d) The N.F.’s and corresponding disc velocities 
are given by the ordinate and abscissa, respectively, 
of these intersections. 


+o 
+— 
| 
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Fic. 3. Natural frequencies and mode of motion for rigid rotat- 
ing disc mounted on a cantilever. 
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(e) The corresponding modes are found from Fig. 
3b. 


A comparison of Eqs. (11) and (13) shows that they 
are both polynomial of the fourth degree in t/a. How- 
ever, Eq. (13) can be solved as a quadratic in 7?/a’. 
Therefore, the labor involved in obtaining all the roots 
of Eq. (13) (for a given value of m) is at the most no 
more than half that needed to obtain all the roots of 
Eq. (11) (for a given value of w). 

If m > A/C, Eq. (13), considered as an quadratic in 
7?/a?, will have one negative root. This negative value 
is discussed because it has no meaning (see Appendix 
I). It arises because of the fact that two of the curves 
in the (w/a)-(r/a) diagram have asymptotes of the 
form 


R(w/a) = t/a (c) 


In Fig. 3a, which is the upper half of the (w/a)-(7/a) 
diagram, the upper curve is seen to have an asymptote 
of the form (c) with k equal to 2. A brief description 
of the behavior of the N.F.’s as w varies will now be 
given. in this description gi is assumed constant. 

Assume that the disc is not rotating. Then the N.F. 
and configuration for the first mode of motion is shown 
on Figs. 3a and 3b by Gy and Mj, respectively. As w 
increases from zero, the moment applied by the disc 
to the shaft decreases and the angle gx decreases. This 
decrease in ga necessitates an increase in the radial 
force, Mr?qio, required to maintain gi constant. There- 
fore, 7 increases. However, as w increases still further, 
go must continue to decrease, and it is obvious that 
go must approach zero asymptotically. This means 
that the radial force, Mr*qi, approaches a finite value 
asymptotically. Therefore, + cannot exceed a given 
maximum limit; for this system the limit is 2a. As w 
decreases from zero, the disc moment causes gz to in- 
crease. The increase in ga necessitates a decrease in 
the radial force to maintain gi constant. Therefore, 
7 decreases. As w approaches minus infinity, go and 
t approach 2(¢.0//) and 0, respectively. 

The N.F. and configuration for the second mode of 
motion, when w equals zero, are given in Figs. 3a and 3b 
by G, and Hp, respectively. An increase in w, from 
zero, decreases the disc moment. (The disc moment is 
less than zero when w equals zero for this mode.) 
Hence, the radial force necessary to maintain gi fixed 
increases. Therefore, 7 increases. As w approaches 
plus infinity, r also approaches the same value. For 
extremely large positive values of w, the disc, in this 
mode, has the same N.F.’s that this disc would have if 
it were rotating with the same angular velocity, w, 
about a fixed pivot at its c.g. with no restoring moment. 
As w decreases from zero, go gradually increases. 
Hence, the radial force decreases and 7 also decreases. 
Since gs must approach zero, 7 must also approach a 
fixed value. Of necessity this value of 7 is the same 


as the upper limit for the first mode of motion. 
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Fic. 4. Angular momentum vectors for a rigid rotating disc 
mounted on a cantilever. 


X 


Fic. 5. Nonrotating rigid disc mounted on a cantilever. 


At the point G; in Fig 3a and H; in Fig. 3b, the.N.F. 
and mode for the first mode of motion are 7 = a and 
q20/quo/l = */2, respectively; the shaft speed is 


A/C (14) 


This means that at the speed given by Eq. (14) the 
system in the first mode behaves exactly as if the disc 
were a point mass. Consider Fig. 4. The plane of 
the paper is assumed to be the OZO’ plane of Fig. 2- 
The angular momenta of the system corresponding 
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Fic. 6a. Natural frequencies for rigid disc mounted on a canti- 


lever. 


| A TT 
— 
| 
| 
a 
Fic. 6b. Natural frequencies for rigid disc mounted on a canti- 
lever. 


to the angular velocjties given above are the marked 
vectors. It is seen that the gyroscopic moment Cwrg2 
and angular inertia moment — Argo are just equal and 
opposite. Therefore, the shaft does not exert any mo- 
ment on the disc for these conditions. 

Consider the system shown in Fig. 5. D is a non- 
rotating rigid disc mounted rigidly on the cantilever 
shaft OO’. Except for the moments of inertia of D, 
this system is exactly the same as the one shown in 
Fig. 2 (i.e., the mass of D’, the length OO’, the flexural 
rigidity of the shaft, etc., have the same values as those 
used previously). Since the disc, D’, is not rotating, 
no value need be assigned to the polar moment of in- 
ertia. The moment of inertia of the disc about a line 
through O’ perpendicular to the plane of the figure is 
assumed to be (A — mC) (where A, m, C have the same 


| 


Fic. 6c. Natural frequencies for rigid disc mounted on a canti- 
lever. 


value as above). It is assumed that the motion of the 
system takes place in the plane of the figure. qo repre- 
sents the maximum linear displacement of the disc 
from the equilibrium position; 2» represents the maxi- 
mum angular displacement of the disc axis from OZ. 
The equation of motion for the system are, after assum- 
ing small vibrations present: 


12B B 
= — — Qo + 
(d) 
B 4B 
— (A —mC)r*q20= — 


If the last of Eqs. (12) is substituted in Eqs. (7) the 
following equations are obtained: 


12B 6B 
—Mr*qwu = — — + — 
6B 4B 
—(A—mC)r*qo= qi — 


It is obvious that Eqs. (15) and (d) are identical. 
Therefore, the rotating system shown in Fig. 2 may be 
replaced by the hypothetical nonrotating system of 
Fig. 5. It is interesting to note that the angular 
moment of inertia of the nonrotating disc becomes 
negative for m > A/C. When this occurs, it is no 
longer possible to use a real disc. 

The above analysis shows that the N.F.’s of the sys- 
tem of Fig. 1 are functions of the three variables w/a, 
a?/8?, and p. Therefore it is not possible to represent 
the N.F.’s on any three-dimensional graph. Instead, 
a family of three-dimensional graphs is necessary. 
Figs. 6a, 6b, and 6c are three graphs of this family. 
The mode of motion, however, is the same for any N.F. 
regardless of the values of w/a, a?/8?, and p. Hence, 
Fig. 3b is complete in itself. 

Fig. 6a plots the N.F.’s as a function of w/a and a?/- 
6°; p = 2in this graph. Figs. 6b and 6c do the same 
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thing except that » = 1 and '/2, respectively. It is 
clear that Fig. 3a is included in Fig. 6a. 

Fig. 6a (as well as Figs. 6b and 6c) shows that as 
a’/B* decreases the system behaves more nearly like 
a point mass. As a?/8? increases, however, the system 
acts more like a disc with fixed c.g. This is as it should 
be because a/ is a measure of the ratio of the radius 
of gyration (corresponding to the polar moment of 
inertia) of the disc to the shaft length. 

A comparison of Figs. 6a, 6b, and 6c (for a fixed but 
arbitrary value of a?/8*) shows that as p decreases 
(from 2) the effect of the angular velocity w also de- 
creases. Ultimately, the angular velocity is of no im- 
portance and may be disregarded. Again this is in ac- 
cordance with physical intuition. 

The common point of the lower set of curves occur- 
ring in each of these figures always has the coordinates 
A/C, 1. As mentioned above, the shaft does not exert 
any moment on the disc at these values of w/a, t/a— 
i.e., it behaves as a point mass for these conditions. 

If a rigid rotating disc is mounted on any other shaft 
arrangement, a new series of graphs must be made. 
This is unfortunate but it cannot be avoided. 

In conclusion, it must be emphasized that the above 
graphs can be easily obtained by replacing the variable 
w/a with m(r/a). This substitution makes it possible 
to solve the characteristic equation as a quadratic in- 
stead of a quartic. The saving of labor entailed in this 
step is obvious. When the number of discs is equal to, 
or greater than, 3, the advantages in using this substi- 
tution are even greater. 


II. N RiGrm Rotatinc Discs ON AN ELAstTic SHAFT 


A. Elastic System Pos;zesses Symmetry 


Fig. 6 shows a system consisting of rigid rotating 
discs D,, Do, ..., D, mounted on a vertical elastic shaft 
S. The shaft is assumed massless, and it may have 
any bearing arrangement. The bearings may be either 
flexible or rigid. When the system is at rest (as shown 
in the figure), each disc c.g. and axis is assumed to coin- 
cide with the straight elastic axis of S. The angular 
velocities we, wi, ..., Ww, of the discs about their axes 
are assumed to be related to the angular velocity of the 
first disc by a constant multiple—i.e., 


o, = Rim 2,4,..., 2n) 


If the discs all rotate with the same angular velocity, 
then they may be considered as rigidly connected to a 
rotating shaft; otherwise each disc can be considered 
as centered on a nonrotating shaft by means of a fric- 
tionless rigid bearing. A turbogenerator set is an illus- 
tration of the former case, while a fan connected to an 
electric motor by means of a belt or gear drive illus- 
trates the latter. Finally, the system is assumed de- 
void of all friction. 

Let C2;, Az; Mo; represent the polar moment of in- 
ertia, diametral moment of inertia, and mass, respec- 
tively, of D;; let cu, C2, ... represent the coefficients 
of elasticity of the assumed symmetrical elastic system. 
Then the equation from which the N.F.’s of the system 
can be derived is 


— C12 C1, Qn 
C29 — (Az — Cokem)r? C2, 2n 
0 (16) 
| | 
| Con, 1 Con, 2n — (Aon — ConRonmyr? | 
where The N.F.’s of the system are best presented in the 
ms a form of a graph similar to Fig. 3a. A graph of this 
(e) type may be constructed as follows: 
w2. = mT (f) (1) Assume a value for m. 


m has the same meaning as in the first section of this 
paper. The mode of motion for any N.F. can be de- 
termined from the linear equations 


(Cn Myr?) Fi + + ..... +61, 2nFon = 0 
(uF; + [G2 — (A2 — Cokom)r?] Fo +..+ C2, nFon = 0 


Can, iF, + Cony + 
+ — [An — Con Ron Fon = O 
(17) 
Appendix I discusses at length the method used in ar- 
riving at Eqs. (16) and (17). 


(2) Find the positive roots, r?, of Eq. (16) using 
the m of (1). 

(3) Take the positive square roots of the values 
found in (2). Denote these by 7, 72, .. 

(4) Using the m of (1), plot the number pairs (rym, 
71), (Tom, T2), ..., (tTg¢mM, Tz) With 7m as abscissa and 7 as 
ordinate. 

(5) Assume a new value for m and begin at (1). 


By choosing suitable values for m between — © and 
+e and following this procedure, 2” curves will be 
obtained in which the abscissa is the angular velocity 
of the first disc and the ordinate is the N.F. 
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Fic.7. mnrigid rotating discs mounted on an elastic shaft. 


The N.F.’s for any given yalue of w, are obtained in 
the following manner. Let w, = wo. Find the inter- 
sections of the curves with the lines w. = wy) and a, = 
—w. Let these intersections be denoted by 71, 72, ..., 
To, and 71’, 72’, ..., Ton’. Then the N.F.’s are 7, 72, ..., 


, , , 


The mode of motion for wo and 7; (say) is found from 
Egs. (17) in the usual manner. However, it is not pos- 
sible to construct a graph similar to Fig. 3b for this 
case. The difficulty is due to the complication involved 
in representing a complete mode as one point on a 


graph. Asa result, the mode is’usually found only for 
a natural frequency that is of particular interest. 

In Section I, the important types of exciting forces 
were classified. Further, a procedure was given there 
for determining when these exciting forces would be in 
resonance with the system of that Section. The only 
modification necessary in that procedure to make it 
applicable here is to relate the exciting forces to the 
angular velocity of D,. If this is done and if a graph 
of the N.F.’s as a function of the angular velocity of D, 
has been obtained by the above steps (1), (2), (3), (4), 
and (5), the procedure described in Section I can be 
used to determine resonance conditions. 

The form of Eqs. (17) suggests that a system with 
nonrotating discs could be devised which would have 
the same N.F.’s as the original. To do this, the only 
change necessary in the system shown in Fig. 7 (which 
has rotating discs) is to alter the value of each 
disc diametral moment of inertia; the altered values 
are 


A, — mk,C, (r = 2, 4, ..., 2m) 


As long as A,/k,C, > m, forr = 2,4, ..., 2n, it is at 
least theoretically possible to construct an equivalent 
mechanical system with nonrotating discs. If A,/k,C, 
< m for any value of r, a mechanical construction is not 
even theoretically possible, because negative or zero 
moments of inertia cannot be represented by real discs. 
However, it is possible to construct an electronic circuit 
that has the equivalent N.F.’s of the original system. 
The negative moments of inertia are analogous to power 
input devices. 

The system with nonrotating discs mentioned above 
can be used to calculate the N.F.’s by either Mykel- 
stad’s tabular procedure* or by Southwell’s Relaxation 
Method.‘ Myklestezd’s method can be applied in the 
usual manner, the negative masses causing no difficulty. 
The Relaxation Method also requires no modification 
beyond avoiding the negative characteristic values. 
This can be accomplished by starting each stage of 


‘relaxation with an assumed configuration that makes 


T; > 0. 


B. Elastic System Unsymmetrical 


The natural frequencies for the system shown in Fig. 
6, with unsymmetrical elastic properties, can be de- 
termined from the roots of 
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The mode of motion for any N.F. may be found from 
Eq. (61) in the usual manner. Appendix II contains 
the derivation of both Eqs. (18) and (61). 

Eq. (18) is a polynomial of 4m degree in r?. It has 
4n real roots, some positive and some negative. As 
before, the negative roots have no meaning and may be 
discarded. 

The method described in IIA for plotting the N.F.’s 
as a function of w,. is applicable for this case. Since it 
can be shown that Eq. (59) has 4” positive real roots 
for any given set of we, a, ..., @2,, this graph will con- 
tain 4n curves. The procedure given in IIA, for de- 
termining the N.F.’s for any given set of we, w,, ..., Wan 
and for finding the disc speeds and N.F.’s at which 
resonance occurs (with different types of exciting 
forces), can also be applied here without modification. 


- The modes of motion cannot be plotted on a graph be- 


cause of their complex character. Therefore, it is rec- 
ommended that the modes be found only for those 
N.F.’s that are of interest. 

The conclusions of Appendix II show that South- 
well’s Relaxation Method can be used to determine 
the N.F.’s and normal modes of the system. As before, 
the negative characteristic values can be avoided by 
always choosing trial configurations that make 


M\(Fy? + Fo?) + Ao(Fiz? + Foo”) + ..... 
2n 
Aon(Fi, + Fa, — m Coker Fir Far >0 (g) 


Myklestad’s tabular procedure can be modified to apply 
to Eqs. (61). However, the rotating discs cannot be 
replaced by nonrotating discs as in the case that pos- 
sessed elastic symmetry. 


APPENDIX I 


Fig. 7 shows n rotating rigid discs mounted on a ver- 
tical elastic shaft, S. Each disc is mounted rigidly on 
the shaft, and it is assumed that each disc c.g. and axis 
coincides with the elastic axis of the shaft when the 
system is at rest. we, ..., are the angular veloci- 
ties of Di, D2, ..., Dn, respectively, about their axis. 
Fig. 8 shows the instantaneous position of the disc, 
D;. The X YZ axes, with origin at O, are fixed in space; 
OZ is vertical and it coincides with S when the system 
is at rest. The X,’Y,’Z,’ axes, with origin, O,’, at the 
c.g. of D;, remain parallel to the X YZ axes during mo- 
tion. 1, 2; — 1 and ge, 2; — , are the x and y coordinates 
of O,’ projected on the XOY plane. qi, 2; and qo, 2; are 
the coordinates of P,’ projected on the X,’0,’'Y,’ plane, 
and they are dimensionless. As in the first part of this 
paper, O,’P,’ is along the axis of D, and its length is one 
unit. @1, 23 — 1 92, 23 — 1, 91, 27, and Qe, 2; are assumed 


small. With regard to the physical characteristics of 
D;, M2; - 1 is its mass, C2; is its polar moment of inertia, 
and Ag; is the remaining principal moment of inertia. 
It is assumed that the shaft, S, and its bearings have 


Z; Pi 


Z 
X 
Y 
Fic. 8. Instantaneous position of the disc of Fig. 7 


equal elastic properties in the XOZ and YOZ planes 
and that it is massless. 

' The kinetic and potential energies for this system 
are: 


2T = Mi(gu? + gar?) + + da”) +... 
+ 2n + Qe, 


— — — (19) 
— Crnwon(Gi, 2nG2, 2n — 91, 2nQ2, 2n) 
and 
2V = ¢ngu’* C1, 2nQi, 2nQ11 
+ + + + C1, one, 
+ Con, 2 + Con, 2nQ1, 
+ Can, 102102, 93 + Con, 2nG2, 2n” 


respectively. The coefficients c,, are the usual coef- 
ficients of elasticity (c,, = ¢s,). T and V are both posi- 
tive definite quadratic forms in the velocities and dis- 
placements, respectively. The form of T and V shows 
that the system has cyclic coordinates whose velocities 
We, W4, Wan- 

The usual procedure at this point is to employ the 
method of “‘ignoration of the coordinates’’® and derive 
the equations of motion in terms of the palpable co- 


° 
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ordinates and the corresponding cyclic momenta. In 
turn, these equations can be used as the basis for a dis- 
cussion of the N.F.’s of the system. However, it will 
be found that at a natural frequency the phase of 
each coordinate is not the same. This means that 
normal coordinates do not exist as in acyclic systems. 
Further, the equations of motion are not a statement of 
the necessary and sufficient conditions that the whirling 
natural frequencies be stationary for slight variations 
in the values of mode. Therefore, this method is not 
followed. Instead, the equations of motion will be 
derived in terms at the palpable coordinates and the 
cyclic velocities. 

The equations of motion, obtained from Eqs. (19) 
and (20) by means of Lagrange’s equations, are 


Migu + Cugu + + .......... + C1, 91, 2» = 0 
Myger + 11921 + +... + C1, 292, 2» = 0 
Cugu + + + +..+ C2, 2» = 0 
+ — + +..+ C2, onG2, n = 0 


(21) 
and 
| coos (guts duge) | = 0 
i 2W2 12922 12922 


Com (G1, 2nG2, 22 — 1, 2n92 | = 
It is obvious that if periodic external forces and mom- 
ents act on the system, the right-hand side of Eqs. (21) 
and (22) would not be zero. In this discussion, how- 
ever, all external forces and moments are excluded. 
The equations of motion given by Eqs. (21) correspond 
to the palpable coordinates and Eqs. (22) to the cyclic 
coordinates. Eqs. (22) show that the cyclic momenta 
of the system are constant. 
By means of the substitutions, 


= M; 17? C12 
C21 (C22 + Cunt — 


Eq. (27) is a polynomial of degree 4n in r. If its 
roots are all real, they determine the N.F.’s of the sys- 
tem. The nature of these roots will now be established. 

Let o be a solution of Eq. (27) and let Fi, Fa, ..., 
Fyn be the corresponding variables that satisfy Eqs. 
(26). Then, 


2n 
+ ( CaF) =0 (29) 


where 


1947 


Q2r—-1 
Qer 


the 4n Eqs. (21) may be reduced to the 2” equations 
Migs + engi + . C1, nGon = O 
+ — + part C2, = O 


Con, 191 + Aandon — + Con, = O 


Qi, 2r + 192, 2r 


It should be noted that this simplification can be car- 
ried out only if the shaft, S, and its bearings have equal 
elastic properties in the XYOZ and YOZ planes. When 
calculating the N.F.’s of the rotating elements in air- 
craft reciprocating engines and gas turbines, this condi- 
tion is usually satisfied. However, there are systems 
that do not satisfy this condition. This method of 
analysis does not apply in those cases. Since Eqs. 
(24) are equivalent to Eqs. (21), the solution of Eqs. 
(24) and (22) follows. 

Eqs. (24) are linear differential equations with con- 
stant coefficients. Therefore, assume solutions of the 


form 
q = Fie" 
(25) 
Substituting Eqs. (25) in Eqs. (24), 
— Mir?) Fi + ....... + C1, = 0 
Cu Fy + (C22 + Cowet — Aor?) Fo +. .+ C2, nF on = 0 (26) 
(Can, + — Fin = 0 
Eliminating the 2n — 1 ratios F,: Fy: .... :F2, from 
Egs. (26), the characteristic equation is 
A(r) = (27) 
where 
C1, 2n 
C2, 2n (28) 
“+ — Angr?)| 
2T0(F) = MiFi? + + ......... + (30) 
2V(F) =enFi? + + ...... + 61, nF iF on 


Both Eqs. (30) and (31) are positive definite quadratic 
forms. The coefficient of the first degree term in Eq. 
(29) may be either positive or negative. Eq. (29) 
shows that, if any of the roots of Eq. (27) are real, they 
may be either positive or negative. Let o’ be another 


root of 
variabl 
its F’s 


2V(F, 


2V0(F,J 
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root of Eq. (27) and F,’, Fo’, ..., Fon’ the corresponding 
variables. Then these quantities, together with o and 
its F’s satisfy 

2n 


2V.(F,F’) + (> — 2T)(F,F’)o? = 0 
2,4 


2n 


2V0(F,F’) + (= o’ — =0 
2,4 
(32) 
In these equations 


27)(F,F’) = + +...+ AonFonFon’ 


2V(F,F’) = tn Fi F;’) C1, onl 
+ C2n, 1F Fo,’ + ee + C2n, Fon’ 

(33) 


If ¢ and o’ are distinct roots of Eq. (27), the following 
relations are valid: 


2n 

oto’ = Cw,F,F,'/2To(F,F’) (34) 
2,4 

oo’ = Vo(F,F’)/To(F,F’) (35) 


Assume that Eq. (27) has an imaginary root, and let 


0, o’ = p— 06 9 
F, = + F,’ = — 0, (r = 1,2, ..., 2m) 
(36) 
Substituting these relations in Eq. (35), 
oo’ =p? + = — Vo(u) — (37) 


To(u) + To(v) 


This is impossible since 7Jo(u) and Vo(u) are positive 
definite. Therefore, Eq. (27) has 4m real roots, and 
these roots are continuous functions of the cyclic veloci- 
ties. This proves that the system shown in Fig. 7 
has 42 N.F.’s for a fixed but arbitrary set of we, ws, ..., 
Won. 

Eq. (29) shows that no roots of Eq. (27) change sign 
as the cyclic velocities vary. However, when all 
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Up to this point nothing has been said about Eqs. 
(22). However, since it has been shown that all the 
roots of Eq. (27) are real, it is clear that Eqs. (22) are 
satisfied by Eqs. (25). Therefore, the cyclic momenta 
of the system are constant. 

The geometric interpretation of this method of find- 
ing the N.F.’s is as follows. Let 


(r = 2,4, ..., 2m) (38) 


k, = 


If the mechanical system under consideration contains 
only geared connections among its rotating elements, 
these ratios are constant. In what follows, these ratios 
will be assumed constant. Using these ratios, we may 
plot the 4n N.F.’s as the ordinate with w: as the abscissa 
and obtain 4” continuous curves. Thus, the above 
process for determining the 4m roots of Eq. (27) is 
equivalent to obtaining the intersection of these 4n 
curves with the vertical line, w. = constant. : 
Examination of Eq. (29) shows that the coefficient 
of the linear term is analogous to the “dissipation 
function” introduced by Lord Rayleigh.’ This ana- 
logy is algebrpic, however, and has no physical signifi- 
n 


cance, for > Cww,F,2| represents an imaginary 


“dissipation function.’’ Since a real dissipation func- 
tion affects the amplitudes of vibration by an exponen- 
tial factor, it is clear that an imaginary dissipation 
function must affect only the periods of vibration of the 
motion. Thus, this concept of imaginary damping 
could be used to prove that the roots of Eq. (27) are 
all real. The next step affects the elimination of the 
linear term in Eq. (29), and introduces the desired hy- 
pothetical acyclic system. 
Let 


(39) 


Substitute these equations 


= MT 
where m is a real number. 
in Eqs. (26). This gives 


(ey F, Cir the wie ote onl’ on = 0 
+ [C2 — (Ae — Fo +..+ C2, on Fon = 0 


Con, iF; + [Con, (Aon | Fon = 


cyclic velocities are zero, 2” of the roots of Eq. (27) are (40) 
positive and 2m are negative. Therefore, 2n of the 
N.F.’s are positive and 2n are negative. The characteristic Eq. (27) now becomes 
— C1 C1, 2n 
2 
C21 C22 (Az Cokom)r C2, 2n =0 (41) 
Con, 1 Con, 2n (Ao, 


This equation is a polynomial of degree 2m in r*. (In 
what follows, Eq. (41) will be considered only as a 
polynomial of degree 2m in 7?.) The transformation 
(39) has thus removed all terms of odd degree from the 


characteristic equation. The nature of the roots, r’, in 
Eq. (41) is the next question. 

Let o? be a solution of Eq. (41); let Fi, Fo, ..., Fan 
be the corresponding set of terms which satisfy Eqs. 
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Using these values, 
(42) 


(40) when o? is substituted for 7?. 
Vo(F) — Ti(F)o? = 0 
where 


27T,(F) = (Ag Crkem) Fe + “oe +- 
(Asn — Conkonm) Fo, (43) 


Now 7;(F) is a function of m; in particular, half of its 
coefficients may be negative or zero. Therefore, T,(F) 
may be either positive, negative, or zero. From this it 
may be concluded that, if Ai(r?) = 0 has real roots for 
7, they may have any values. 

Assume that o? is another root of Eq. (41) distinct 
from o?; let F,’, Fo’, ..., Fon’ be the corresponding 
These quantities, in conjunction o?, Fi, Fo, ..., 
Fy, and Eqs. (40), permit the formation of 


Vi(F,F’) — Ti(F,F’)o? = 
Vi(F,F’) — Ti(F,F’)o’? = 0 


Hence, if o? and o”? are distinct, it may be concluded 
that 


(44) 


Vo(F,F’) 0, Ti(F,F’) =0 (45) 


In this equation 


= + (A, — +... 
+ (Aan — Conkonm) FonFon’ (46) 


These [Eqs. (45)] are the usual orthogonal relations ob- 
tained for conservative acyclic systems. They depend 
only on the symmetry of 7; and Vo and not on the 
signs of the coefficients of these two forms. 

If Eq. (41) has imaginary roots, let 


o* = p + 4, 
F, = pp + try, \ (ry = 1, 2, ..., 2n) 
(47) 


Substituting these relations in Eqs. (45), the following 
equations result: 


= p — 46, 
FP,’ Mr + 1, 


Vo(u) + = + 
Ti(u) + Ti(v) = 0 


The first equation of Eqs. (48) is impossible, since Vo(F) 
is a positive definite quadratic form in the variables F 


(48) 


(or »). This proves that the characteristic equation, 
Eq. (41), has no imaginary roots. 
Consider the function 


U = T,(F)/Vo(F) (49) 


Vo(F) is a positive definite quadratic form and 7;(F) is 
a quadratic form whose coefficients are either positive, 
negative, or zero. Since Vo(F) does not vanish (except 
when F; = F; = ... = Fo, = 0), it is obvious that, for 
at least two sets of Fi, Fo, ..., Fn, U has a stationary 
value. If the method of the Calculus is used to search 
for one of these values, it is seen that Eqs. (40) are the 
required conditions and that the reciprocal of one of 
the roots of Eq. (41) is the value in question. Thus 
U has at least one real stationary value; it may be 


either positive or negative. Let Fi, Fo, ..., Fo, be the 
corresponding set of variables. If the search for sta- 
tionary values is continued,’ it can be shown that 2n 
stationary values exist and that these are the reciprocals 
of the 2n real roots of Eq. (41). Therefore, the station- 
ary property of the roots of Eq. (41) has been estab- 
lished. 

By means of a linear transformation of the variables, 
the two forms 7;(F) and Vo(F) may be reduced to 


271(6) = + 
2V.(8) => + 


This transformation is always possible as long as Vo(F) 
is positive definite. Hence, normal forms of 7;(F) and 
Vo(F) exist. In addition, the following relations are 
valid: 


vr = 1/o,2 (r = 1,2, 2n) (51) 


where 017, 02”, ..., 2” are the 2 roots of Eq. (41). 

A theorem of Algebra*® states that the number of 
positive y,*’s is equal to the number of positive coef- 
ficients in Eq. (43). Thus, by determining the values of 


A, — C,k,m (7 = 2,4, ..., 2n) (52) 


we can determine the number of positive y,”’s. Inas- 
much as Eq. (43) is a function of m only, it is apparent 
that there can be no less than m and at most 2n positive 

It was shown above that the system under consider- 
ation has 4m N.F.’s for a given set of w’s. Further, if 
these N.F.’s are plotted as ordinates with w as abscissa, 
4n continuous curves result; 2m are on the positive side 
of the horizontal and 2n are on the negative side. The 
introduction of Eq. (39) thus means that we are deter- 
mining the intersections of these 4n curves with the line 

= mr (53) 
Thus, the presence of negative roots of 7? in Eq. (41) 
or of negative 7,” has no physical meaning and merely 
shows that some of the N.F. curves have asymptotes 
of the form of Eq. (39).* Therefore, the square roots 
of the positive roots, r?, of Eq. (41) are the N.F.’s of the 
system. 


* A simple algebraic proof of this can be given as follows: 
Let o; be a root of Eq. (27) for a given set of w’s. 

Then find m, ko, ky ..., Ren, so that m = w2/o;, and ky = w,/ws 
(r = 2,4, ..., 2m). Substitute these relations and r? = o;? in 
Eq. (41). A comparison shows that Eqs. (27) (with r = o;) and 
(41) are identical. Hence, any one of the 4” roots of Eq. (27) is 
the square root of one of the positive roots of Eq. (41). Let 
a2” be one of the positive roots of Eq. (41) for a given set of m, 
Ba, 

Find we, ws, ..., @2n SO that we = oom, and w, = kywe (r = 4, 6, 
..., 2n). Substitute these relations and + = az in Eq. (27). 
It is seen that Eqs. (27) and (41) are identical. Hence, the 
square root of any of the positive roots of Eq. (+1) is also a root 
of Eq. (27). Therefore, there is a one to one correspondence: 
between the 4n roots of Eq. (27) and the square roots of the 
positive roots of Eq. (41). The negative roots of Eq. (41) 
(thought of as a polynomial of degree 2n in 7*) have no meaning. 
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APPENDIX II 


If the system shown in Fig. 6 (and discussed in Ap- 
pendix I) does not have symmetrical elastic properties 
the details of Appendix I are changed but the procedure 
used there remains the same. Therefore, in this sec- 
tion only the details changed by the assumption of 
unsymmetrical elasticity will be pointed out. 

With different elastic properties in the XOZ and 
YOZ planes (see Fig. 7), the elastic potential energy 
becomes 


2V’ =cngu® + cugugu +...... + C1, 2n 
+ + + ..... + C1, 2n’Qe, 
Con, 191191, 2n « Con, 2n Qi, on? 


The coefficients of elasticity satisfy 


Corrs Cre = Cor 


Cn = 


The altered equations of motion are 


Migu + Cugu + + ..... + C1, 291,21 = 0 
Migr + + +..... + C1,2n'Ge, 2 = 0 
+ + Cowege2 + Coogie + ........ =0) (55) 
+ — + +........ = 


This equation is a polynomial of 8 degree in r. 
It can be shown that Eq. (58) has 8m real roots, 4n 
positive and 4m negative. Hence, the unsymmetrical 
system has 4 positive and 4m negative N.F.’s for any 
given set Of we, ws, ..., Won. If all are related to we 
by a constant multiple, it is possible to plot a graph 
similar to Fig. 3a with w. as abscissa and N. F. as ordi- 
nate. In this way a graph is obtained with 4m curves 
above the zero ordinate and 4n below. 
By means of the equations 


w, = kymr 


This equation will be considered as a polynomial of 
degree 4n in 7?. 
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Eqs. (22) remain the same. 


Assume Eqs. (55) have solutions of the form 


qu = Fu sin zt 
qu = Foy cos rt 
qi2 = Fis sin tt 
Q22 = Fe cos tt 


Substituting Eqs. (56) in Eqs. (55), 


> Myr?) Fa + C12’ Foo = 
Cu Fy > (Coe Aor?) Fie Foo +. eee 


oooo 


(56) 


(57) 


It is easily seen that Eqs. (56) satisfy Eqs. (22) regard- 


less of the values of Fo, ..., Fe, on, 
Eliminating the 4n — 1 ratios Fy:Fu: ... 
from Eqs. (57), the characteristic equation is 


A’(r) = 0 


eee 


Substituting these equations in Eqs. (57), 


Myr’) Fu + + ..+ a, 
— Mir?) Fe + G12’ Foo + ..+ 61, on’ Fo, = 0 
iu (C22 Agr’) Cokomr? Foe + 
+ Fi, = 0 
Foy (C29’ Aor?) Fre Cokymr? Fy. 
+ C2, Fe, m = 0 


a simplification can be made in the above analysis 


Fa, 


(58) 


It may be shown, by the same procedure that was 


used in Appendix I, (1) that Eq. (62) has 4m real roots; 


17 
where 
M;,r? 0 C12 0 
0 éu’ — 0 C12" 
0 C2, — — Cowet (59) 
0 Ca’ — Cxwet C22’ — Agr? 
(61) 
7 fs (7 = 2, 4, ...» 2n) (60) The characteristic equation, for this case, is 
0 Cy = M. 17? 0 C12’ 
Ai’ (7?) = C21 0 C22 — — Ckemr? = 0 (62) 
0 C2 Cokomr? Cas! Aor? tere 
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(2) that some of these roots of (r?) may be negative, 
and (3) that the roots of Eq. (62) are the stationary 
values of the function 

2n 


oT" 
2V0'(F) 


In Eq. (63) 


2T0'(F) = Mi(Fir? + Fair®) + Fis? + Feo”) + ..... 
+ Mon — — 1 + Fe, 17) + Aan(Fi, an? + Fe, 2n”) 
(h) 


2V0'(F) = euFu? + t+. . + C1, ut 
C11’ Foy? + C12! Foo Fait. . +61, 2n’ Fe, 
Cm, Fuki, + Con, nF i, + 
Con, 1' Fai Fo, + Con, 2n’ Fo, on? 


However, it is not possible, in this case, to find a simple 
method for determining the numbers of negative roots 


which Eq. (62) possesses for a given set of ko, hi, ..., 
Ron, m. As in Appendix I, the negative roots of Eq. 
(62) have no meaning and the positive roots determine 


the N.F.’s. 
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Analytic Definition of a Retractable 
Landing Gear Axis of Rotation 


ROY A. LIMING* 


North American Aviation, Inc. 


SUMMARY 


This paper presents an analytic solution for the problem of 
determining a skewed axis of rotation for a retractable landing 
gear whose motion involves a single articulation. 

Two simple linear slope equations, developed to define the 
axis, provide the foundation required for the study of the kine- 
matics system. 

The application of this analytic technique reduces the solution 
of a given problem to a few minutes of mechanical manipulation 
of a standard calculating machine. 

While there are other procedures, both graphical and mathe- 
matical, available for obtaining a solution to the problem, the 
author believes that the method suggested represents the simplest 
in principle and in practical application to design, as evidenced by 
the fact that the discussion in this paper represents an abstract 
of the design procedure used in the North American Aviation 
Engineering Department with respect to landing gear de- 
velopment in the preliminary as well as the production design 
stage. 

Two examples are worked out in complete detail to show the 
high-speed inspection character of the analytic approach. 


INTRODUCTION 


¢ THE DEVELOPMENT of a retractable landing gear the 
designer is usually confronted with the initial 
problem of establishing an axis of rotation that will in 
one articulation move the wheel from a given extended 
to a given retracted position. The discussion that 
follows is concerned with the general case, in which the 
plane of symmetry of the wheel is displaced, thereby 
introducing the familiar problem of the skewed line in 
space geometry. 

The determination of this axis of rotation, a perennial 
problem child in aircraft design circles, has been the 
subject of a wide variety of attacks, both graphical and 
mathematical. A careful survey of such attempts re- 
veals that the mathematical approach has been selected 
almost invariably as an incidental or subordinate corol- 
lary to some particular graphical technique. The 
special criteria that should properly govern such 
analysis have been largely ignored, thereby placing the 
mathematical technique at a disadvantage. 

The analytic solution presented in this paper results 
in the development of two simple linear slope equations 
defining the axis of rotation. The extreme simplicity 
of procedure and high speed which characterize the 
establishment of these slope-intercept equations actu- 
ally revolutionize design procedure to the extent that 
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all graphical layout intended to define the axis is ren- 
dered superfluous. 


DESIGN PROCEDURE 


This analytic method is employed directly to define 
quickly and accurately the axis of rotation as soon as 
the designer has selected wheel positions, extended and 
retracted, which represent a reasonably satisfactory 
design configuration. 


The usual design approach is to first establish the ex- 
tended position of the wheel in the side view. The 
stowed position is then achieved by selecting the loca- 
tion of the center of the wheel in such a manner as to 
reach a satisfactory compromise with adjoining struc- 
ture and the airplane contour. This latter condition 
may require tilting the wheel axle through one or two 
angles of inclination, e.g., either in the front or profile 
view or in both. The tread, variable within specified 
limits, is selected to fit this general pattern. 


The coordinates of the center of the wheel, extended 
and retracted (see A, and A, of Fig. 1), are thus given. 
Any second point in the rotating system may be em- 
ployed as a basis for calculation. For maximum con- 
venience, this point may be located a distance of unity 
along the centerline wheel axle from the wheel center, 
either inboard or outboard. The coordinates of this 
point (see B, and B, of Fig. 1) may be established 
immediately in terms of the axle position. The axis of 
rotation may then be precisely determined with an 
average expenditure of 10 min. time, using Chart 1 and 
a standard calculating machine. 


Should the axis of rotation thus determined fall in 
such a position within the airplane as to prevent prac- 
tical design of a trunnion to support the gear, it becomes 
necessary to relocate the wheel, usually in the extended 
position, with respect to the tread and to calculate a 
new axis position. If the designer, because of the 
stowage problem, finds his degree of freedom largely 
restricted to the tread variable, he may request the 
simultaneous application of the maximum inboard, 
mean, and maximum outboard positions of the wheel 
extended. The trend of the rotation axis, subject to 
the tread variable, may thus be established at once. 
If no adjustment is necessary, the axis stands precisely 
defined by the slope-intercept equations obtained by the 
use of Chart 1. (See Appendix.) 
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Landing Gear Axis of Rotation. 


CuHartT 1 


Analytic Calculation Chart 


Yi 
Center\4p: 33.3 30 —24.9 
Xe Zs 
wheel An 0 0 
Point x Y3 
centerline ) BE: 32.3 30 —24.9 
w Xs Ys Zs 
axle (Br: —0.190067 0.254101 0.948318 
Direction numbers of AEAR: (X1 — X2):(¥i — ¥2):(Zi — Z2) 
33.3 30 —24.9 
Direction numbers of BEBR: (Xs — Xs):(¥a — Ya):(Za — Zs) = a: be 
32.490067 29.745899 —25.848318 
A M2 | qa Bi Bz 
—51.746321 —0.672060 | 0.306040 2,628.9 2,562.3 40.113 0.856 
2A = —103.492642 | 


Plan view equation: 
Side view equation: 


X = MY + Bi:X = —0.672060Y + 40.113 
Z = + = 0.306040Y + 0.856 


CRITERIA CONSIDERED 


The following criteria have been established to 
govern the analytic approach to this problem: 

(1) Recognition of analytic techniques as the natural 
mathematical tool for the solution of point-line-plane 
relationships involved in motion analysis in space. 

(2) Employment of a systematic analytic attack to 
effect a completely simplified solution of the general 
case. 

(3) Elimination of all intermediate steps and formu- 
las between the original basic data and the actual proc- 
ess of calculating the required final data. 

(4) Reduction of actual calculation to the most 
efficient possible mechanical routine. 


© WHEEL AXLE 
(RETRACTED) 


AXIS OF 


WHEEL AXLE 
(EXTENDED) 


Fic. 1. Perspective view of relevant areas of the wing show- 
ing the basic geometric definition of a retractable landing gear 
axis of rotation. The accompanying discussion analytically 
reduces the kinematics dimensional specification of this design 
problem to a simple algebraic routine. 


(5) Recognition of the superiority of analytic tech- 
niques in the subsequent development of the complete 
dimensional specification of the landing gear structure 
and its coordination with the cambered surface of the 


aircraft. 


GEOMETRIC DEFINITION OF THE PROBLEM 


Fig. 1 includes relevant plan, profile, and elevation 
areas of a typical wing combined in convenient perspec- 
tive form. To establish the geometry of a retractable 
undercarriage, it is necessary first to define, within the 
available stowage space, a rotation axis that satisfies 
the two given positions (extended and retracted) of the 
wheel. 

The wheel is defined geometrically by a straight line 
representing a segment of its own axis of rotation (the 
centerline of the axle). The geometric problem con- 
sists, then, of bringing this straight line from its ex- 
tended position in space to its retracted position within 
the cambered surface of the wing. 

The axis of rotation that will permit the achievement 
of this motion may be defined as a skew line in space, 
determined by the intersection of two planes normal to, 
and bisecting, the two lines determined by two positions 
of any two given points in the rotating system. The 
basic principle involved is shown in Fig. 2. For the 
geometric solution of this problem the reader is referred 
to any standard text in descriptive geometry. 


ANALYTIC SOLUTION 


The analytic solution is analogous to the geometric 
definition of the axis. In Fig. 1, the center of the wheel 
(A) has been chosen as one point, and an arbitrary point 
(B) conveniently located a distance of unity along the 
centerline of the axle from point (A) has been chosen as 
the second of two points defining the position of the 
wheel, the subscripts (E) and (R) denoting the extended 
and retracted positions, respectively, of the two points. 
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Fic. 2. In this yiew, the rotation axis is seen as the axis of 
acone of projection. The loci of all points of the rotating system 
are circular arcs in this view and yield the same angular displace- 
ment for any given angle of retraction. This view also provides 
an edge-view of the perpendicular bisector planes of Fig. 1, and 
supplies from elementary geometry the proof of the method used 
to define the axis of rotation. 


The origin of the coordinate system employed may be 
taken at any convenient point—e.g., at the conventional 
fuselage reference system origin, or at some translated 
point, such as the intersection of a perpendicular 
through the wheel center with the fuselage plane of 
symmetry. In the practical example solved in the 
Appendix, this translated origin is established at the 
center of the wheel retracted (see point A, of Fig. 3). 
The justification for this choice will be obvious to the 
reader. 

In the selected system of coordinates the given points 
are defined as follows: 

Center of wheel 


Ag: Yi, Zi) 
Ap: (Xo, Y2, Z2) 
Point on centerline wheel axle 
Bg: (Xs, Y3, Z3) 
Br: (X4, 


The next step is the analytic description of planes per- 
pendicular to lines AzAr and BzBy at the respective 
midpoints of the two lines. 

Coordinates of the respective midpoints of the two 
lines are as follows: 


Xit+X. Vit Y2 


Li 
ine AgAR 9 2 


Direction numbers of the two lines are as follows: 
Line (X1 — Y2):(Z; — Ze) = 
Line BzBp: (Xs X4):(¥3 Y,):(Z3 Zs) = 


Then, from elementary analytic geometry, the equa- 
tion of the plane perpendicular to line AgA, at its mid- 


RETRACTABLE LANDING GEAR AXIS OF ROTATION 


point (plane I of Fig. 1) may be written directly as 


x(x - %) +(¥%,— x 


and of the plane perpendicular to line ByB, at its mid- 
point (plane = of Fig. 1) as 


— (x + + (Ys — Ys) X 


(2) 


Eqs. (1) and (2), when transposed, rearranged, and 
simplified, become 


(Xi — X2)(Xi + X2) 4 


aX = 


— + (Z; — Z2)(Z; + Z2) (3) 
2 2 
(Ys — VYs)(Vs + Va) + (Zs — Z.)(Zs + (4) 


2 
or, more simply, 
+42 = + Vi? + — 
Y,* — Z;*) (5) 


+ + = + Va? + — — 
Y? — 2,7) (6) 


Eqs. (5) and (6), treated as a simultaneous system, 
define the line of intersection of the two planes, i.e., 
the required axis of rotation. 

A simpler set of equations is obtained by referring to 
the direction number concept of analytic geometry for 
more convenient slope-intercept description of the axis. 
Since the axis appears most closely related to the Y-axis 
in a position sense, a direction number pattern of the 
type M,:1:M, applies,’ yielding slope-intercept equa- 
tions as follows: 


X=MiY+ 
Z= M:Y + 


To achieve this pattern, Z and X are eliminated, 
respectively, from the simultaneous system .consisting 
of Eqs. (5) and (6), whence 


bet, — dice ap— ag . 
Y 7 
X + oA (7) 


dab, — ng — (8) 
A 2A 


in which A = — 
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FRONT VIEW 


Fic. 3.@For illustration convenience, the 1-in. dimension locating the point Bz (or Br) is applied to exaggerated scale. 
designer may desire to apply a slight angle of ‘‘toe” to the wheel, thereby modifying the X and Z coordinates of point Bz. 
pose the wheel is rotated counterclockwise (in the front view) about its center (Az) through an angle a. 
X, — cos a, Z; = Z; + sin a, used in the chart, will impose this condition upon the kinematics pattern. 


Xe + Ve + (Xe + Ve + 


ll 


and 
gq = X3? + Y3? + — (XP + Ye + 


CONCLUSIONS 


Eqs. (7) and (8) constitute the only working equa- 
tions* required to define absolutely and precisely the 
axis of rotation. These equations supply immediately 
the tangent of the angle that the axis of rotation makes 
with the Y-axis in both the plan and side view, as well 
as the coordinates (B,, 0, B.) of one point on the line. 
Additional points may be quickly obtained by direct 
substitution in these linear expressions. 

Chart I indicates the small number and complete 
simplicity of the calculation steps required to secure 
complete analytic definition of the axis of rotation. 

In the example illustrated in Chart 1, based upon an 
actual design problem, only seven separate calculation 
steps are required, on a calculator featuring cumulative 
addition and subtraction of products, to secure the 
required slopes and intercepts. Only three of these 
steps are incidental; four out of the seven represent the 
desired final slope-intercept data defining the axis of 
rotation. 

The mechanical routine of Chart 1 consumes an 
average of but 10 min. time, a fact easily checked with 
the data given in the example. Maximum precision 
may be observed with no additional time required, 
since it is convenient to use a calculator to the limit of 


Ac (Be) 


SIDE VIEW 


The 
Sup- 
Then the values X; = 


its keyboard, simply by recording more decimal places 
during the calculation process. 

When the four basic points defining the position of 
the wheel have been set by design, the analytic pro- 
cedure described above will greatly expedite the pre- 
liminary design stage of development, as well as dictate 
the dimensional specification of the final design, to the 
complete exclusion of cumbersome, inaccurate graphical 
procedures. 


APPENDIX 


The high degree of practicability of the two working 
equations, Eqs. (7) and (8), is fully demonstrated by 
the solution given below for the kinematics problem 
suggested in Fig. 3. 

The data shown represent the actual basis for suc- 
cessfully stowing the main wheel of a tricycle gear 
within the cambered surface of a thin profile, high- 
speed, laminar-flow wing, the “rigging” of which in- 
volved an angle of extreme forward sweep. 

Solution: From Fig. 3, the coordinates of points Ag, 
Ap, and Bg, a point 1 in. inboard along the center- 
line axle from Ag, may be tabulated directly in the 
chart. 

The coordinates of Bz may be obtained by the follow- 
ing analytic inspection method: 

Step 1. Establish from Fig. 3 direction numbers of 
the centerline axle in the retracted position. If a pro 
jection of unity is assumed on the Z-axis, these numbers 
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— tan 11° 20’:tan 15°:1 


or 


— 0.200425 : 0.267949: 1 


Step 2. Calculate the square root of the sum of the 
squares of these three numbers. Divide this value 
into each of the numbers. The resulting quotients are 


— 0.190067 :0.254101 :0.948318 


These values are the direction cosines of the center- 
line axle retracted; since the point By is 1 in. from the 
origin (Ag), they are likewise the coordinates (X4, V4, 
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Z,) of the point Bz and are tabulated directly in the 
chart. 

Step 3. Perform the mechanical calculation routine 
indicated in the chart. (Note that, with Az as the 
origin, the expression X;? + Y,* + Z,? equals 0 in 
solving for p; the expression X,? + Y,? + Z,? equals 1 
in solving for g.) 

REFERENCES 

1 Liming, Roy A., Practical Analytic Geometry with Applications 
to Aircraft; The Macmillan Company, New York, 1944. 

? Liming, Roy A., Airfoil Development and Structure Within the 
Cambered Surface, Aero Digest, December, 1945; March, April, 
June, 1946. 


Dear Sir: 

In the papers “‘Jet Propulsion and Its Application to High- 
Speed Aircraft,’’ by D. J. Keirn and D. R. Shoults, JouRNAL oF 
THE AERONAUTICAL SCIENCES, August, 1946, and “‘Some Ele- 
mentary Details of Ramming Intake and Propulsive Jets for 
Airplanes,” by S. A. Moss and W. R. Foote, JOURNAL OF THE 
AERONAUTICAL SCIENCES, March, 1946, two different interpre- 
tations are given for propulsive efficiency although the equations 
are the same. 

In the first paper it is stated that maximum output occurs when 
% = 2V, at which point the propulsive efficiency is computed to 
be 66.7 per cent. Maximum efficiency is stated as occurring 
when v; = V where output is zero. 

In the second paper one may find a formula for the coefficient 
of propulsive efficiency which shows that the best speed for an 
airplane with jet propulsion is infinite. 

The question arises whether a definition that gives maximum 
efficiency when the output is zero or one that permits efficiencies 
to become greater than 100 per cent is logical from a physical 
standpoint. 

Probably in no other branch of engineering are there so many 
definitions for one and the same coefficient as there are in this 
case. It is proposed to present here the derivation of a formula 
for the propulsive efficiency which is logical from a physical stand- 
point and is not subject to the criticisms that have been made 
above. . 

The general expression for an efficiency is: 


(1) 


In our case we may state more in detail that 7 is the ratio of the 
energy imparted to the body to that of the maximum possible 
amount of energy available at our disposal in the system con- 
sidered. 

The mass of fuel utilized per second, W;/g, possesses, just before 
the combustion process starts, a certain heat content that enables 
it to propel the jet unit at a momentarily constant forward veloc- 
ity, V. In the thermal jet, sufficient air pressure must be built 
up at the entrance to the combustion chamber to force the exhaust 
jet to discharge to the rear. Assume that the mass of air has a 


= output/input 


negligible velocity with respect to the propelled unit as it enters 
the combustion chamber. The heat content, or chemical poten- 
tial energy, will then be used to accelerate the mass of air and fuel 
from zero velocity with respect to the propelled unit to the jet 
Thus 


exhaust velocity, 


Letters to the Editor 


Ly = [(Wa + Wy)/2g]0;* (2) 


where L, is the kinetic energy input to the jet from the chemical 
potential energy and is dependent on the thermal efficiency condi- 
tions whereby the chemical potential energy is converted into 
kinetic energy of the jet. In our discussion we assume the coef- 
ficient of thermal efficiency to be equal to 1. 

But the chemical potential energy is not the only energy that 
the fuel has. Because of the forward velocity of the propelled 
unit, V, the fuel also contains kinetic energy equal to: ‘ 


Lz = (W;/2g) V? (3) 


Thus the maximum energy that can be imparted to the jet unit 
as input under the most favorable conditions is: 


Input = LZ; + Lz (4) 


Since output can be expressed as input minus the losses, con- 
sideration will now be given to the inherent losses in a jet pro- 
pulsion unit. 

Mention has been made of the necessity of bringing the air to 
rest at the entrance to the combustion chamber. The work done 
in compressing this air is chargeable to the thermal efficiency. 
The kinetic energy of the entering air relative to the propelled 
body 

L; = (Wa/2g) V? (5) 


is a loss of kinetic energy as far as the propulsive jet is concerned. 
It is true that it is converted to potential energy in the form of 
pressure and tends to improve the thermal efficiency of the unit, 
but from the standpoint of the total balance of kinetic energy it is 
a loss. 

Similarly, the kinetic energy of the jet with respect to the pro- 
pelled unit differs from the kinetic energy with respect to the 
ground assumed as the reference system. This reduces the ef- 
fective output of the jet by the amount: 


Ly = [((Wa + Wy)/2g](0; — V)? (6) 
Ly is usually called the kinetic energy of the exhaust gases. Ex- 
pressing the output of the jet as the input minus the losses gives: 


Output = La + ~ ~ | ( 
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The reactional force of the outflowing gases is equal to: 


ZL W. + W,; Wa 
pa (tM, y 
V g g 
If the instantaneous coefficient of the propulsive efficiency 7, of a 


jet propulsion unit is to be defined as output divided by input, we 
have: 


(8) 


L 2[(We + — 


7, = const., vj = const., we obtain from 


(9) 


For = const., 
(dn,p/dV) = 0, 
= (We + W;)/(Wa + (Wa? + WaW; + W,?)*/?] (10) 


for 
= [Wa + (Wa? + + + Wy) (11) 


These formulas may be transformed for two extreme cases: 
(a) Pure rocket, i.e., without taking the air from outside: 


W.=0 F= P= vjV (12) 
g g 
np = 2V0;/(V? + (13) 


W, 
Npmaz. = 1 for V = Pmaz. = = (Li + Ly) (14) 


where P denotes the energy output. 
(b) Neglect the mass of the fuel: 


W=0, F= P= - as) 


np = — V)V/o;* (16) 


"pmaz. 


Wefv;\? 
= for = 2V or V = = AG} = 
(17) 


The result is that the coefficient of the propulsive efficiency in 
the case (b) cannot be equal to 1. The coefficient 7, represents 
well the real conditions. For V = 0, P = np = O because the 
body is at rest and one cannot talk about a “propulsive” effi- 
ciency. It has a maximum value of 0.5 for V=0;/2 (v; is assumed 
to be constant) when P achieves its maximum value. It is equal 
to zero for V = vj—i.e., where the momentum at the entrance 


and at the exit is equal. 
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Similarly, the value of the power P, cannot ever utilize the en- 
tire kinetic energy of the jet (W./2g)v;?, since maximum P is 
equal to (Wa/g)(v;/2)?. 

Let us point out that the above definition of the instantaneous 
coefficient of the propulsive efficiency represents the useful per- 
centage of the total power without ambiguity. In the papers 
cited above the coefficient 7, is given identically by the ratio: 


? V2) + 


From this identical equation in one paper it is concluded that 
maximum efficiency occurs when v; = V where output is zero and 
that for the maximum output v; = 2V and np = 66.7 per cent. 
(Actually, the maximum output is equal to 


(Wa/g)(v;/2)? = (Li/2) 


with W; = 0 in Eq. (1), (2), and mp is 50 per cent and not 66.7 
per cent.) In the other paper, after some transformations, it is 
concluded that the best speed for an airplane with jet propulsion 
is infinite and that the propulsive efficiency is equal to 1 when 
vj = V and the thrust and the net power used to drive the air- 
plane are both equal to zero. One may ask what this “propul- 
sive’ efficiency represents if it is not equal to zero and both 
the propulsive force and propulsive power are equal to zero, 
Thus, for the same coefficient that should represent to a certain 
extent the useful power, we have interpretations that represent 
anything else but the useful power. 

Undoubtedly, the best speed of an airplane occurs when the 
optimum value of the useful power may be achieved. A simple 
differentiation of the expression for the useful power (performed 
by Mr. Keirn and Mr. Shoults) shows that the maximum occurs 
for v; = 2V. Hence, this is the best speed. The ambiguity in 
both cited papers has its origin in the fact that all the authors 
consider the losses only in the form of the energy of jet discharge 
gases. (The above derivation of the instantaneous coefficient 
of the propulsive efficiency appeared in the paper ‘Optimum 
Height of Flight of Winged Rocket- or Jet-Bomb,” by M. Z. 
Krzywoblocki, Journal of Applied Physics, Vol. 16, No. 5, pp. 
275-279, May, 1945, and in other magazines.) 

Z. KRZYWOBLOCKI 
R. W. McCroy 
Associate Professors 
Department of Aeronautical Engineering 
College of Engineering 
University of Illinoi 
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Dear Sir: 

I should like to make a few comments on the paper ‘‘On the 
Steady Flow of a Gas Through a Tube with Heat Exchange or 
Chemical Reaction” by Paul Chambré and Chia-Chiao Lin which 
appeared in the Journal of October, 1946. 

(1) The interesting conclusion is reached that, in the steady, 
frictionless, nonadiabatic flow of a compressible fluid through a 
channel of constant cross section, the temperature of the fluid is 
reduced by the addition of heat to the fluid for Mach Numbers 
between 1/+/y and 1. This same result was obtained by the 
writer some years ago, but it was not published. It was also 
obtained and published! by Dr. Bruce L. Hicks of the Cleveland 
laboratory of the N.A.C.A. 

A slight extension of the analysis and a figure can considerably 
clarify the physical picture of the phenomenon. Conservation of 
energy is expressed by Eq. (5) of the reference paper. By the use 
of Eqs. (1), (4), and (14) and the relation 


dp/p = (dp/p) + (dT/T) 


Eq. (5) can be rewritten as 


dQ = + (cp — + udu = dU +dW+ dE 
where U = internal energy per unit mass, W = flow work done 
by unit mass, E = kinetic energy per unit mass, and other sym- 
bols are as in the reference paper. 

For dn = 0, Eq. (22) becomes 

aT/T = [(1 — yM*)/(1 — M*)]|(dQ/H) 

or, by use of H = cyT, 
dH/dQ = (1 — yM*)/(1 — M?) 


It can then easily be shown that 
dE/dQ = (y — 1)M*/(1 — (a) 
dU/dQ = (1/y)((1 — yM*)/(1 — (b) 


dW/dQ = [(y — 1)/y] — — (c) 


Heat energy that is added to the fluid is therefore divided, in 
the proportions given by Eqs. (a) to (c), among a change in kinetic 
energy of directed motion along the pipe, a change in kinetic 
energy of random motion, and a change in flow work or pumping 
work. The proportional division of the heat energy among these 

(Continued on page 63) 
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SUMMARY 


This paper contains a simplified presentation of the exact and 
approximate solutions of supersonic flow over two-dimensional 
bodies. The data are all presented as a variation with the free 
stream Mach Number, thereby enhancing their usefulness. An 
exact solution is given for the maximum flow deflection that 
¢an be obtained by means of an oblique shock wave at a given 
Mach Number. The third and higher terms are presented for 
the infinite series approximating oblique shock flow. The third 
term obtained is not the same as that originally given by Buse- 
mann in 1935 at the Fifth Volta Congress. The new third term 
is analyzed to show that under certain conditions an oblique 
shock flow may be considered isentropic even though the third 
order of the flow deflection or pressure discontinuity is not 
negligible. A reliable graphical method is presented whereby 
the pressure distribution may be quickly and accurately obtained 
for a two-dimensional body in a supersonic flow field. 


INTRODUCTION 


. jew PERIOD has passed when supersonic flow was 
merely of theoretical interest to aeronautical 
engineers. Supersonic wind tunnels are used exten- 
sively and designs for supersonic aircraft are being 
planned. 

This paper is a summary of the problems of two- 
dimensional supersonic flow, an extremely important 
field, since it can be shown that in some cases the flow 
over an entire wing in supersonic flight may be con- 
sidered two dimensional. 

All of the computations are based on a y value of 
seven-fifths or 1.400. The author’s reasons for this 
selection have been previously stated! and are in agree- 
ment with the recent evaluations by Keenan and 
Kaye* based on existing spectroscopic data for the 
constituents of dry air. The thermodynamic deriva- 
tions are based on the perfect gas relation p/p = RT 
and the assumption that the heat transfer is negligible 
and the processes are adiabatic. The additional as- 
sumptions made in the analysis of shock flow are that 
the conditions are everywhere uniform ahead of the 
shock wave, that the equation of conservation of 
momentum may be directly applied to a stream tube 
passing through the discontinuity forming the shock 
wave, and that the loss in kinetic energy through the 
shock wave appears as heat energy that may be directly 
measured as an increase in enthalpy. The flow through 


Presented at the Aerodynamics Session of the Fourteenth 
Annual Meeting, I.A.S., New York, January 29-31, 1946. 
* Section Head of Flight Research Engineering. 
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the shock wave is then an irreversible adiabatic process 
leading to an increase in entropy. 

The standard exact and approximate solutions of 
two-dimensional supersonic flow problems are extended 
to larger flow deflections and higher Mach Numbers. 
The results are made more ‘useful by being presented 
as a variation with the free stream Mach Number. 
A graphical method is presented which will quickly 
and accurately determine the pressure distribution 
about a two-dimensional body in a supersonic flow 
field. This method immediately shows whether or 
not a known solution exists for any portion of the 
body as long as the surface considered is inclined at 
less than 15° to the free stream velocity vector. In 
addition, the method quickly determines a suitable 
profile shape for any given physically possible pressure 
distribution. 


SYMBOLS 


a = local rate of pressure propagation and local velocity 
of sound = vpb/p, (ft. per sec.) 
Cy = specific heat at constant pressure = 6,007.1, 
(ft. per sec.)?/( °F.) 
C, = specific heat at constant volume or constant density 
= 4,290.8, (ft. per sec.)?/(°F.) 
h = enthalpy = + (ar = 
P/} jo 
ho + | pe 
y = 
sec.)? 
H = stagnation or total absolute pressure, (Ib. per sq.ft.) 
M = Mach Number = V/a 
Mimin. = minimum Mach Number for existence of oblique 


shock 

p = absolute static pressure, (Ib. per sq.ft.) 

= pressure coefficient = (p — po)/'/2p0Vo* 

R = gas constant = (Cp — C,) = (53.345)(32.174) = 
1,716.3, (ft. per sec.)?/( °F.) 


if 2 
S = entropy = 5+ G oe ( b/p ): (ft. per sec.) 


Po/po” (°F.) 
absolute Rankine temperature = (°f + 459.69) = 
(°F.) 
velocity component normal to shock, (ft. per sec.) 
velocity component parallel to shock, (ft. per sec.) 
magnitude of resultant velocity vector, (ft. per sec.) 
isentropic gas constant = C,/C, = 7/5 
angle of stream tube deflection or angle of velocity 
vector behind shock (J) relative to velocity 
vector upstream of shock (V;), rad. 
= maximum possible angle of stream tube deflection 
for a given free stream Mach Number, rad. 


er 
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° 
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timum 
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= angle made by shock wave relative to velocity vector 
upstream of shock (V;), rad. 
p = mass density = (Ib. sec.?/ft.*) 


Subscripts 


0 = free stream conditions at an infinite distance upstream 

1 = conditions in uniform parallel flow upstream of shock 
wave which are equivalent to the free stream condi- 
tions at infinity 

2 = conditions in uniform parallel flow behind or down- 
stream of shock wave 

H = conditions at stagnation point where the velocity is zero 


ISENTROPIC FLOW 


In the first section the equations for isentropic or 
reversible adiabatic flow will be developed. Although 
this is an irrotational frictionless flow with no head 
loss, it is necessary for a study of oblique shock waves 
since the flow is isentropic on either side of the ex- 
tremely thin discontinuity sheet defining the shock 
wave. It will also be shown that even for angles of 
flow deflection of 15° the flow is approximately isen- 
tropic across the shock discontinuity. 

The energy equation, neglecting gravity and other 
extraneous effects, is 


V?/2 + h = const. (1) 


This expression of constant energy is true for any heat 
insulated flow; that is, it holds for any adiabatic flow 
whether it is reversible (isentropic) or irreversible 
(friction losses) as long as there is no heat transfer to 
the fluid and the heat produced by any friction losses 


remains in the fluid. 
For a reversible adiabatic or isentropic flow the 


enthalpy may be simply expressed as 
hy — Io = So'dp/p (2) 


Then the energy equation for isentropic flow be- 
comes 


(Vo2/2) + S°dp/p = (Vi2?/2) + S'dp/p (3) 


Because of the fact that the entropy is constant, the 
pressure and density are explicitly related by 


p/p” = const. (4) 


and, since p/pT is constant for a perfect gas, the 
physical states are related as follows: 


Pr Pi Ti a;? ©) 


Consequently, if only the relations at the end points 
of an isentropic flow process are desired, Eqs. (3) and 
(4) may be combined directly in the form 


Vo? y po _ 
2 po 2 y—1p1 


Or, more generally, it may be stated that throughout 
the entire field of isentropic flow 
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2 y—1lp 


For zero velocity the stagnation or total pressure H, 
is obtained, and therefore 


2 y—-lp y—-lp 2 


23: & (7) 
y¥—1 po 


where H is constant as long as the flow is isentropic. 


This may be written as 
4 a’? ay? Vo? 


wat. guy 


Therefore, the Mach Number may be written as 


a y¥—1\a? 


or from Eq. (5) 


(9) 


Fig. 1 shows the variation of p/H with Mach Number 
as obtained from Eq. (9). 

The following useful relations for the static pressure 
may be obtained from Eqs. (5) and (7) 


y-1 V2 

¥/(y¥-1) 

Po PH To Tx 


(10) 


Eqs. (9) and (10) apply to either subsonic or super- 
sonic flow as long as there are no shock waves in the 
region being considered. For subsonic flow the 
determination of the local velocity or pressure is 
difficult since the flow depends not only on the condi- 
tions the stream tubes have passed but also upon the 
conditions everywhere ahead of the stream tubes. 
This is'so because the partial differential equation of 
subsonic flow is of the elliptic type, which has no real 
characteristic values. Consequently, the boundary 
values everywhere in the flow field (including infinity) 
to some extent influence all the flow. The problem 
is comparatively simpler for supersonic flow since, 
in this case, the partial differential equation governing 
the flow is of the hyperbolic type, which does have real 
characteristics.*4 These characteristics are the Mach 
lines, which are defined explicitly by the local Mach 
Number (J) and are straight lines inclined backward § 
at the angle sin~'(1/M) to the local velocity vector. 
The velocity component normal to the characteristic 
or Mach line is always equal to the local. speed of 
sound (a). The physical significance of these charac- 
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teristics is that all disturbances in the flow field are 
propagated along them; that is, the characteristics 
are the envelopes produced by sound waves moving 
at the local speed of sound. A discontinuity, generally 
called a limit line, develops whenever these charac- 
teristics intersect and most shock waves may be de- 
fined as the envelope formed by these intersecting 
characteristics. In general, nothing occurring down- 
stream of a characteristic can directly influence the 
flow field upstream of the characteristic. 

The remarkable advantage resulting from the 
existence of these characteristics for isentropic super- 
sonic flow computation is that as long as no discon- 
tinuities or shock waves are formed, the pressure and 
velocity changes depend only upon the angle the 
stream tube is deflected (@,) and the free stream Mach 
Number (Mo). Consequently, any given solution for 
the pressure change due to a stream tube deflection 


‘at a given free stream Mach Number may be applied 


to any other isentropic supersonic flow as long as the 
initial and final angles of flow are known. The simplest 
solution for isentropic flow is the particular case of free 
expansion around a sharp corner. Then the rate of 
expansion is the maximum physically possible and the 
characteristics or Mach lines become radii originating 
from the sharp corner; consequently, the velocity 
component normal to any of these radii is always equal 
to the local speed of sound, and along a given radius 
the velocity, pressure, and density are constant.5 The 
equation for irrotational supersonic flow expansion 
(—6,) starting from a free stream Mach Number of 
unity is then 


= 


tan— +! tan 
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[oro 


The values obtained are given in Table 1. 

The numerical evaluation for any other supersonic 
isentropic flow may now be obtained by reading @, as 
the difference in Table 1 from the 6, value corresponding 
to the given free stream Mach Number. A compres- 
sion would correspond to a positive flow deflection 
(+@,) while an expansion corresponds to a negative 
flow deflection (—6,). In applying the solution to a 
given surface, however, there are three conditions that 
can invalidate the results. Since the isentropic ex- 
pansion given in Table 1 is the most rapid physically 
possible, the flow would leave a given surface and 
continue exactly as in Table 1 as a free expansion if 
the magnitude of the negative surface deflection were 
greater than the corresponding value given in Table 1. 
The other two limitations are for supersonic compres- 


TABLE 1 


Values for Isentropic Expansion Around a Sharp Corner 
Two-Dimensional Flow 


y = 7/5 
P/H —6,° M 
0.5283 0 1 
0.5 0.437 1.046 
0.45 1.986 1.132 
0.4 4.139 1.223 
0.35 6.789 1.322 
0.3 9.935 1.433 
0.2724 11.87 1.500 
0.25 13.65 1.559 
0.2026 17.81 1.700 
0.2 18.06 1.709 
0.15 23.49 1.897 
0.12 27.44 2.040 
0.1 30.53 2.157 
0.08 34.28 2.300 
0.06 38.75 2.484 
0.05 41.45 2.601 
0.04 44.62 2.746 
0.03 48.50 2.935 
0.02722 49.75 3.000 
0.02 53.61 3.208 
0.015 57.00 3.406 
0.01 61.49 3.693 
0.005 68.45 4.210 
0.001 81.77 5. 566 
0.0005 86.51 6.234 
2.356 X 10-5 102.32 10 
1.278 X 10-8 121.23 30 
2.790 X 127.59 100 
0 130.45 


sion. The flow must not become subsonic from too 
great a compression. The rate of surface deflection 
must not be large enough to make the characteristics 
intersect so near to the surface that the resultant 
shock wave would affect the surface pressures.® 

When the values in Table 1 are plotted in the manner 
shown in Fig. 1, all three of the previous objections 
are automatically considered for any given problem. 
For example, the @, scale is also shown in Fig. 1 for 
any isentropic flow starting from a free stream Mach 
Number of 1.5. Then if the negative surface de- 
flection for free expansion is greater than 118.6° the 
Mach Number becomes infinite, showing that the 
flow has left the surface and has continued as a free 
expansion. If the positive surface deflection is greater 
than 11.87° then the Mach Number becomes less than 
unity, showing that the solution is invalid. 

So far the solution is exact for any flow expansion and 
for extremely small flow compressions. The solution 
is also exact for the particular case of large flow com- 
pressions where the characteristics or Mach lines do 
not intersect. However, it will be shown later that 
even when the Mach lines do intersect and form shock 
waves, Fig. 1 still gives a good approximation to the 
flow pressures for positive surface changes of less than 
15° since the flow is still approximately isentropic. 
If shock waves are then admitted to the solution 
another difficulty arises. If +06, is greater than 
Osmaz. (Appendix I), there is a standing shock wave 
upstream of the surface producing mixed subsonic and 
supersonic flow behind the shock wave and over the 
surface. However, in Appendix I it is shown that for 
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any +8,° less than 15° the use of Fig. 1 will auto- 
matically yield a subsonic Mach Number if the standing 
wave exists in front of the surface. Consequently, 
Fig. 1 not only provides an extremely rapid solution 
for supersonic flow problems now solved by existing 
methods but it is also remarkably reliable. 

The pressure change for a given flow deflection may 


be obtained from Fig. 1 by 


Pr _ 
/H 


ISENTROPIC 
COMPRESSION 


+ 
©, + 
Ce, 


PRESSURE COEFFICIENT 
6 


610s 


3 
C,0,+ 
_ISENTROPIC 
5 EXPANSION 
Fic. 3. 


and the pressure coefficient P may be computed as 


Vi? (y/2)M;? Pi 


or for isentropic flow 


1 1 + [(y — 1)/2] 
+ [7 - 1) 


Fig. 2 shows the variation of the pressure coefficient 
with a given stream tube deflection starting from a 
free stream Mach Number of 1.5. Fig. 3 shows the 
effect of the initial free stream Mach Number on the 
pressure coefficient produced by a stream tube de- 
flection of 5°. 

The approximation to the isentropic solution as found 
by Busemann* yields 


P = CA, + + C02 + ..... (12) 


where 


= C= [(y + 1)/2] Mit — 2M)? + 2 


(My? — 1)?’ (Mi? — 1)? 
1 (y + 1) 
C3 = 4 Mi ies 
+5) ape 
(4y* — 287? + 1097? + 2307 + 105) oars +4] 
48(y + 1) 3 


and the pressure coefficients obtained from these terms 
are also shown in Figs. 2 and 3. 

The first term C,0, is simply the tangent to the 
isentropic solution and is the linear approximation. 
The addition of the second term C,6,? takes into} ac- 
count the curvature of the isentropic solution / and 
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yields the parabolic approximation. The addition of 
the third term C;0,* considers the rate of change of 
curvature of the isentropic solution and for M > 1.3 
gives a good approximation even for large —@, values 
(isentropic expansion) and gives a good approximation 
for +6, (isentropic compression) up to about 5°. 
These approximations are extremely useful for ana- 
lytical solutions. For small surface deflections, general 
conclusions for supersonic airfoils have been obtained 
by formally integrating the expression obtained from 
the first two terms over the airfoil surface.*—® 


OBLIQUE SHOCK WAVES 


The equations for two-dimensional supersonic oblique 
and normal shock flow will now be developed. These 
relations were discovered by Rankine, Hugoniot, and 
Meyer, and although they are available in many 
standard references*—> they will be presented here in 
a form more suitable for actual numerical computations. 

Fig. 4 shows the oblique shock wave discontinuity 
with an abrupt change in pressure and velocity. The 
conditions of uniform parallel flow prevail on each side 
of the discontinuity sheet or shock wave. Although 
the previous characteristic solution for isentropic 
supersonic flow could exist for any positive or negative 
6,, the oblique shock wave exists only for a finite 
positive 0,; that is, only for an abrupt flow compression 
accompanied by an increase in entropy.® 

The momentum equation applied to the flow normal 
to the discontinuity yields 


pis? + Pi = + po (13) 


The continuity equation for the flow normal to the 
discontinuity is 


pilly = polls (14) 
The momentum equation for the flow tangential to the 
discontinuity is 
= (pete) v2 
and comparison with Eq. (14) .shows that necessarily 


OBLIQUE SHOCK FLOW 


/ 


Fie. 4, 


(15) 
The energy equation, Eq. (1), becomes 


or 
(u?/2) + = (ua?/2) + he (17) 
In addition, for any perfect gas 


Therefore, Eqs. (16) and (17) become 


2 
2 — 1 pa 
+ — 
2 — 1p 


— + = (19) 
2 2 y—lp 


Eq. (18) shows that the stagnation temperature (Ty) 
must be the same on both sides of the discontinuity. 
It is interesting to note that Eq. (18) is identical in 
form to Eq. (7). However, it will be shown that the 
following inequalities are associated with Eq. (18) 


> H, < 
PL 


while for Eq. (7) these quantities must both be equal 
in order that Eq. (5) may be used. This means that 
for isentropic flow, the pressure or density is deter- 
mined explicitly from Eqs. (5) and (7), while some addi- 
tional relation between the pressure and density is 
required for the adiabatic flow process before Eq. (18) 
can provide a complete solution. 

For example, if we take an adiabatic flow with a 
head loss but with the additional restriction that the 
local static pressure (2) has returned to the original 
free stream static pressure (po) then 


po _ 
po To a? 
Also, Eq. (18) may be written identical to Eq. (8) as 
Vo? Ae? ag? V2 az? 
2 y¥-i y¥-1 2 y¥-1 


Consequently, for = po and Hy 


2 
ey 
2 Vo? 


ao” p2 
1+ — 1)/2]M,? 
{ [(y )/2] 


1 + [(y — 1)/2]M;? 
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flow with 
pr Po Fh 


may be written as 


p. a,” pe 


- 
on 


{ i+ Ky 

Of course, the latter relations for isentropic flow all 
reduce to unity when p2 = pp. The comparison simply 
serves to emphasize that Eqs. (3) and (5) must be used 
only for isentropic flow. 


For adiabatic shock flow with a head loss and entropy 
increase, the pressure and density relation may be 
obtained directly from Eqs. (13), (14), and (19) by 
eliminating and 


[Cv +1)/(y — +1 _ +1 

+ /(y — 1)) + (b2/b1) + 

br _ [r+ — 1 6(p2/p1) — 1 

pr + 1)/(y — — — 
(20) 

The isentropic flow equation, Eq. (4), is compared with 

the shock flow equation, Eq. (20), in Fig. 5. 

From Fig. 4 and Eqs. (14) and (15) 


Therefore 


(1 sin 6, 1) 
p2 sin 0, cos — 
sin 0, (22) 
COS sin (A. — 85) 
Also from Fig. 4 
um = V,sin 0, = aM, sin 6, (23) 
= V2 sin (0. — = sin (0, — 


The pressure discontinuity may now be evaluated by 
combining Eqs. (13) and (14) in the following manner: 


p2 
Then Eqs. (22) and (23) may be used to eliminate 1 
and p: from Eq. (24) 
sin 6, sin 6, pi Vi? 
cos (A. — 6,) 1 + cot 0, cot 6, 

(25) 

The velocity discontinuity may also be obtained by 
combining Eqs. (20) and (24) 


(24) 


fi = 


From Eq. (10) an analogous expression for iSentropic 
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21 


(26) 
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Then, introducing Eq. (23) 
M;? sin? 


M.? sin? (0, — 0,) = -— = 23 
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Another useful relation for M, in terms of only @, and 
6, may be obtained by writing Eq. (26) as 

Yor 
pi 2 pi 


= 


and the desired equation is obtained by substituting 
Eqs. (23) and (25) 


y+ 1 sin 6, sin @, (28) 


sin? 6. — 
cos (6,, — 


M? 2 


Expressions for the stagnation pressure (7) on either 
side of the discontinuity may be obtained from Eqs. 


(9) and (27) 
(y-1)/ ) 
(4)" = yy 
pr 2 
4y sin? 0, 1 y¥+1 
Hp (y-1)/y | (29) 
2 
4y sin? (6. — pe y¥+1 


Although the discontinuity in entropy across the 
shock wave is associated with an irreversible flow 
process, its magnitude may be evaluated by integrating 
any reversible process between the two states since 
entropy is a thermodynamic property that is inde- 
pendent of the process and is only a function of the 
initial and final physical states. 


30 
(80) 


For any isolated system the entropy must either remain 
constant or increase and the entropy can remain con- 
stant only for a reversible adiabatic, therefore 


(p2/p2”) > 
and from Eqs. (20), (27), and (29) 
po > py > pr, M2 < Mi, < 
Consequently, the shock wave defined by these equa- 


tions is an irreversible adiabatic compression resulting 
in an increase in entropy. 


NorMAL SHOCK WAVES 


The preceding equations also apply to normal shock 
waves when 6, = 7/2 and @, = 0. In order to apply 
Eqs. (13), (14), (16), (18), and (20) direetly to a normal 


shock, it is necessary only to note that 4% = V and 
v = 0; consequently, the relations for a normal shock 


are 
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y—1L\p 


and 
y — 1L\p2 2y hr 27 
or in terms of M, 

_ yMi? — [(7 — 1)/2] 


pi (y + 1)/2 
Pr _ (y + 1)/2 
+ [(y-—1)/2] Ve 
1+ [ty — 1)/2)M? _ (Tx/T;) 
— [(y — 1)/2] [Cv + 1)/2] 
= {1 + [(y — 1)/2) 


pe 4y + 1 


DIscussION OF SHOCK WAVE EQuaTIONS 


If the oblique shock wave relations for a constant 
free stream Mach Number (MM) and a varying flow 
deflection (6,) are desired, then the equations and 
graphs introduced by Taylor and Maccoll® give a 
convenient solution. However, the following pres- 
entation is more useful, especially if the oblique shock 
wave relations are required for a given angle of flow 
deflection as would normally be required for a given 
wedge or airfoil profile. 


Taking y as 7/5 the equations for oblique shock flow 
become 


6 — (p2/p1) 
[(6/7) (p2/p1)] (1/7) M2? 


sin? 6,, 


[(6/7)/(b2/p1)] + (1/7) 
sin? (0, — 


Hy _ + th _ 


H, pili + (2/5) pi 


The curves presented in Figs. 6 to 10 were obtained 
for constant values of @, by assuming values of 6, be- 
tween 0 and 7/2. Also included in the figures are the 
limiting cases of the oblique shock; that is, the normal 
shock wave (6, = 1/2, 6, = 0) 


po 1/6(7M,? 1), 6(p2/p1) +1 


Mg +1) Ha _ bs + 
7 pil 6+ 


and the Mach wave (6, = sin—(1/M;), 6, = 0) 


po tan 6, p2 6(p2/p1) — 1 
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= po/p. = M2/M, = = 1 
The normal shock and the Mach wave lines form an 


envelope for the oblique shock waves as @, approaches 


zero as a limit. Figs. 6 to 10 are equivalent to the 
Busemann shock polars** but are more suitable for 
aumerical computation. 
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It is important to note in Figs. 6 to 10 that there is 
a minimum free stream Mach Number (Mimin.) below 
which an oblique shock wave solution does not exist 
for a given @,, and when this minimum Mach Number 
is exceeded there are two possible oblique shock wave 
solutions for each M;. The solution for Mim. for 
the existence of an oblique shock with a given stream 
deflection is obviously identical to the solution of the 
maximum flow deflection (@smer.) possible at a given 
free stream Mach Number. The exact solution of this 
problem is given in Appendix I and Table 2. The 
numerical values from Table 2 are also included in 
Figs. 6 to 10 and form a line dividing the oblique shock 
flow solution into two parts. The solid line portion, 
although it has a smaller entropy increase, is still the 
more dynamically stable oblique shock wave’ and is 
the one always found in practice.® 

Fig. 9 shows that M2, the Mach Number behind 
the normal shock wave, is always subsonic, and the 
same is true for the flow near Mimi. for an oblique 
shock. However, for the oblique shock M; is only 
slightly less than unity, and as Mimm. is exceeded the 


TABLE 2 
Values of the Minimum Mach Number and Maximum Flow 
eflection 
Two-Dimensional Flow 
y = 7/5 
mes, Minin M2 Pi 
0 1 1 1 1 90 
0.5 1.046 0.9855 1.074 1.053 0.9999 80.28 
1 1.075 0.9775 1.122 1.086 0.9998 77.93 
2 1.122 0.9660 1.205 1.142 0.9994 75.12 
2.5 1.143 0.9614 1.243 1.168 0.9990 74.12 
5 1.240 0.9438 1.482 1.291 0.9952 70.82 
10 1.421 0.9251 1.843 1.537 0.9717 67.45 
15 1.614 0.9188 2.360 1.813 0.9385 65.73 
20 1.84 0.9201 3.073 2.142 0.8695 64.90 
25 2.125 0.9281 4.133 2.546 0.7581 .62 
27 2.265 0.9328 4.722 2.736 0.6989 64.63 
30 2.519 0.9408 5.895 3.058 0.5908 64.80 
35 3.144 0.9566 9.370 3.723 0.3705 65.40 
40 4.446 0.9764 19.179 4.610 0.1306 66.33 
45 14.10 0.9977 198.2 5.829 0.00943 67.62 
45.584 1 6 9 67.792 
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abrupt decrease in Mach Number becomes less and the 
flow remains supersonic throughout. Figs. 6 to 10 
all show that once Mimi. has been exceeded the rela- 
tions for the oblique shock wave tend to form a family 
of curves parallel to the simple isentropic Mach wave 
relation, especially for small @,. It should be possible 
to determine simple approximate relations for oblique 
shock waves; for example, a family of parallel straight 
lines would approximate the relation between M, and 
M, (Fig. 9) for a rather wide range of M; and 4,. 

At present there is no known general flow solution 
downstream of the shock wave if the flow near the 
surface has become subsonic. This is so because, as 
previously mentioned, the resulting subsonic flow by 
itself is rather difficult to solve since there are no real 
characteristics, while the actual mixed field of subsonic 
and supersonic flow is even more complicated. 

For M;, less than Mimim. there is a standing shock 
wave upstream that is entirely detached from the two- 
dimensional body. For a region near the nose of the 
body the flow is subsonic with an ordinary stagnation 
point somewhere on the nose of the body. No solution 
has yet been presented for the pressures over the body 
in this type of flow but Crocco" has found a relation 
for the amount of rotation in each stream tube behind 
the shock wave, and the author, in a paper not yet 
published, has obtained an approximate solution for 
the distance from the standing shock wave to the 
stagnation point on the body. Of course, a body 
having a finite thickness at the nose or having a semi- 
wedge angle greater than 45.58° will always have a 
detached standing shock wave at ordinary Mach 
Numbers. Only a true wedge having @, less than 
Omer. Will allow an ideal attached oblique shock wave 
to exist for a finite free stream Mach Number. 

If M; is only slightly greater than Mimim. the shock 
wave will be attached to the nose of the wedge but the 
flow behind the shock wave will still be subsonic over 
a portion of the wedge and no general solution can be 
obtained. The shock wave will not be a plane dis- 
continuity sheet but will be curved near the wedge 
and gradually will become a plane parallel to the free 
stream Mach lines at some distance from the body. 
The flow through the outer portion of the shock wave 
will pass through unaltered and of course be super- 
sonic as well as isentropic. 

However even if the flow is supersonic everywhere 
behind the oblique shock, a true plane oblique shock 
can exist only for a finite distance from a finite body. 
Only an infinite wedge can have an oblique shock wave 
producing uniform parallel flow everywhere behind 
the shock wave. The infinite wedge and its infinite 


oblique shock wave must of necessity also have an 
infinite drag since the entropy discontinuity, although 
it may be small, still acts over an infinite surface. A 
Mach wave can extend to infinity since it is an infinitesi- 
mal disturbance which has no discontinuity; however, 
a shock wave always has a finite discontinuity. A 
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finite body, of course, cannot have an infinite drag, so 
the oblique shock of a finite body cannot extend to 
infinity. This may also be found to be true from the 
physical consideration of the Mach lines formed by the 
changes in the body’s surface behind the nose wedge. 
These Mach lines interfere with the oblique shock at 
some finite distance from the body. The Mach lines 
form an envelope which causes the oblique shock wave 
to curve and continually reduces the pressure and 
resulting entropy discontinuity until at some finite 
distance from the body the oblique shock has been 
transformed into a Mach wave and the flow is then 
unaltered in passing through it. 

Since the stagnation enthalpy (fy) is constant 
through a shock wave, Eq. (16), Crocco!! has shown 
that the flow rotation in any stream tube behind a 
shock wave is proportional to the static pressure (2). 
Consequently, the flow immediately behind a shock 
wave is irrotational (velocity potential exists) only 
if the flow is uniform and parallel. Only an infinite 
wedge satisfies this condition exactly. However, for 
all practical purposes the oblique shock existing at 
the nose of any body having a sharp leading edge 
whose @, is less than Osmgz, satisfies these conditions. 
The oblique shock will remain plane far enough from 
the body that the flow may be considered everywhere 
irrotational behind the shock if the flow is everywhere 
supersonic. Therefore, a velocity potential exists 
between shock waves and it is valid to use the isentropic 
solution of Fig. 1 for any flow expansion. 


APPROXIMATE SOLUTION OF OBLIQUE SHOCK FLow 


Fig. 10 shows that the head loss is almost negligible 
for flow deflections less than 5° and is fairly small even 
for larger flow deflections at some Mach Numbers. 
This is an indication that the flow is approximately 
isentropic for small flow deflections. That this is so 
was first proved by Busemann’® who found that 


— Si = + 0(6,*) 


However, Figs. 5, 7, and 10 indicate that even a fairly 
large 0, may possibly produce an oblique shock across 
which the flow is approximately isentropic if the pres- 
sure or density change is small. In this case the two 
curves in Fig. 5 can be shown to be identical to the 
second order because the first and second derivatives 
of Eqs. (6) and (20) have the same limit as M, ap- 
proaches unity (p: — pi, p2 — pi). For example, Eq. 
(26) can be written as 


\p2 — Pi P2 \p2 — Pi 


or 


Then, when M, approaches unity as a limit the above 
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equations for an oblique shock becomes identical to the 
isentropic equations 


dp/dp = y(b/p) = a? 


The curves in Fig. 5 diverge with a large pressure or 
density increase because the third and higher deriva- 
tives at an M, of unity are not identical. It is seen 
that if p2/p: < 1.60 or po/p; < 1.95, then the oblique 


shock 


1)'+ 11) (2- 1) = 


+ 1) | 


4 
1 
) 20 
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‘due to the shock will be less than 3 per cent. 


The preceding conclusions may also be justified by combining Eqs. (20) and (30) to show that for an oblique 
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shock flow Eq. (20) is represented within 1 per cent 
by the isentropic flow Eq. (6). Figs. 7 and 8 show that 
for Mach Numbers less than 2.71 a flow deflection of 
10° will be within this range, while for smaller flow 
deflections that Mach Number range within which the 
flow is approximately isentropic is greatly increased. 
Fig. 10 shows that for this range the total pressure loss 


1) - 


3 — 6 cos? 
12 sin? cos? 2 sin 6, cos 6, (31) 


so that if the third order of density, pressure, or flow deflection increment is negligible the oblique shock flow is 
Then if the conditions of negligible entropy change are satisfied, Fig. 1 or Table | 
This is a great simplification since Fig. 1 then provides 


approximately isentropic. 
may be used for the oblique shock flow computation. 
a solution for all supersonic flow problems in this range. 


where be supersonic and that the minimum Mach Number must be exceeded before an oblique shock can exist 
are automatically satisfied since in Appendix I it is shown that Fig. 1 fails to provide a solution whenever these 


conditions exist for any +4, less than 15°. 
For normal shock 


therefore 


for a given positive flow deflection (+94,) 


[(b2/p1) — 1] = [2y/(y + — 1) 


Sz Si = — 1)/3(y + — 1)? + — 


Consequently, even a normal shock is approximately isentropic for Mach Numbers near unity. Figs. 7 and 8 
show that the normal shock can be approximated by isentropic compression within 1 per cent for Mach Numbers 
less than 1.35. Fig. 10 shows a 3.03 per cent total pressure loss through the normal shock at M,; = 1.35. 

In Appendix II an infinite series expansion is found for the pressure increase through an oblique shock wave 


by = 2/(My? — 1), by = {[(y + 1)/2] Mit — + 2}/(My? — 1)? 


In addition, the conditions that the flow must every- 


(32) 


1 + 1? (—3y? + 127 + 7) + 1) 
Ms + 4 — 2M, + = 
P 16 ‘ 12 fo 
1 (3y3 — 3y? — 7y — 1) (3y* — 3y? + 337 + 7) (37? — 9y — 2) 
10 M$ 
(Sy + 2) 2 ] 
3 


The expressions for }; to bs are derived in Appendix II; the values of b; to b, for various Mach Numbers are givel 


in Table 3 and presented in Fig. 11. It is seen that C, and C, of Eq. (12) are identical to b, and & of 
therefore the variation of C; with Mach Number is also shown if 


Eq. (32); C3, however, is not identical to bs, 
Fig. 11. . 
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The first two terms of Eq. (12), or Eq. (32), give 
the two approximations for oblique shock and super- 
sonic isentropic flow which have been used in practice.® 
The first term (0,0, or C,6,) gives the linear or Ackeret® 
approximation while the addition of the second term 
or gives the parabolic or Busemann® ap- 
proximation. The additional term C;0,° giving the 
third approximation to isentropic flow is taken as 
given by Busemann® but the author believes that the 
term (C; — D)0@,* giving the third approximation to 
the oblique shock flow as derived by Busemann® is 
incorrect and the expression derived in Appendix II, 
b0,3 is the correct third order term. 

Fig. 11 shows that for Mach Numbers less than 1.3 
the values for 5, b2, bs; and C; tend to approach infinity 
rapidly as M, approaches unity; consequently, the 
power series representation by Eq§. (12) or (32) is not 
good for Mach Numbers less than 1.3. 

In Fig. 11 it is seen that b; and C; are of approxi- 
mately the same magnitude and that for Mach Num- 
bers greater than 1.82, }; is less than C;. This shows 
that, contrary to Busemann,® the oblique shock wave, 
at least for Small flow deflections, produces a smaller 
pressure increase than does an isentropic compression 
for all free stream Mach Numbers greater than 1.82. 
In addition, the b, term becomes negative for all Mach 
Numbers greater than 2.4, indicating that for larger 
flow deflections at high Mach Numbers the oblique 
shock produces a much smaller pressure increase than 
does an isentropic compression for the same flow de- 
flection. This result should be expected since the 
increase in entropy due to the oblique shock is at the 
expense of the kinetic energy which is converted into 
heat energy (an increase in enthalpy). Consequently, 
for a given speed the pressure change should be less 
as the entropy increases. However, Fig. 11 shows 
that b; is slightly greater than C; in the Mach Number 
range of 1.34 to 1.82. Therefore, for a given free 
Stream Mach Number in this range and for a given 
small positive flow deflection, the oblique shock pro- 


duces a slightly greater pressure increase than does an Large M 
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isentropic compression. This apparently anomalous 
condition can exist because at a given free stream 
Mach Number the pressure increase in Eq. (32) de- 
pends only on the flow deflection. Figs. 7 and 8 show 
that for this Mach Number range the pressure or 
density discontinuity is actually a minimum for any 
given small flow deflection. That is; for a given @, 
the resulting p2/p1 or p2/p; is a minimum and is closest 
to that for a Mach wave; consequently, from Eq. 
(31) the entropy discontinuity is actually a minimum 
in this Mach Number range. In addition, the data 
for Fig. 10 definitely showed that for any given small 
6, the value of H;/H; was nearer to unity for Mach 
Numbers in the range of 1.34 to 1.82. 

It is rather interesting to digress for a moment and 
note that for an airfoil in subsonic flight the local Mach 
Number range in which the compressibility shock! 
occurs is generally 1.4 to 1.5. This is the Mach 
Number range wherein 5; is approximately 9 per cent 
more positive than C;. Consequently, in this local 
supersonic Mach Number range for a given pressure 
increase the required positive flow deflection is slightly 
less for an oblique shock than for an isentropic com- 
pression. Now the positive flow deflection on the rear- 
ward portion of an airfoil section could be produced 
only by a boundary layer separation or thickening. 
Therefore, if the boundary layer can easily leave the 
surface, the required compression for a return to sub- 
sonic flow could occur by means of a relatively large 
boundary layer separation accompanied by an isen- 
tropic compression and Mach Number decrease outside 
the boundary layer” (the boundary layer flow is, of 
course, subsonic and nonisentropic) until finally the 
boundary layer will no longer thicken. The final flow 
compression to a subsonic Mach Number can then be 
accomplished by a normal shock wave. However, if 
the boundary layer cannot easily separate from the 
surface, then the flow compression and resultant sub- 
sonic Mach Number could be obtained directly by an 


TABLE 3 
Values of Coefficients for Power Series Approximations 
y = 7/5 
M by = Q be = Ce bs Cs bg 
1.1 4.364 30.32 544.4 646.3 1299. 
1.2 3.015 8.307 653.22 59.05 516.3 
1.3 2.408 4.300 14.53 14.94 77.49 
1.4 2.048 2.919 6.128 5.885 22.51 
1.5 1.789 2.288 3.331 3.059 8.965 
1.6 1.601 1.950 2.153 1.961 4.293 
1.8 1.336 1.618 1.280 1.262 1.339 
2.0 1.155 1.467 1.019 1.140 0.4711 
2.2 1.021 1.386 0.9394 1.174 0.1274 
2.4 0.9167 1.337 0.9356 1.263 —0.0468 
2.6 0.8333 1.306 . 9654 1.374 —0.13819 
2.8 0.7647 1.284 1.013 1.494 —0.1540 
3 0.7071 1.269 1.069 1.619 —0.2248 
4 0.5164 1.232 1.405 2.252 —0.2928 
5 0.4082 1.219 1.764 2.879 —0.3221 
6 0.3381 1.212 2.127 3.498 —0.3340 
8 0.2520 1.207 2.852 4.720 -—0.3470 
10 0.2010 1.204 3.578 5.937 —0.3471 
100 0.02000 1.200 37.28 59.99 —0.350 
2/M 1.2 0.36M 0.6M -—0.352 
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oblique shock accompanied by a relatively small sepa- 
ration of the boundary layer. Figs. 6 and 10 show that 
for a given pressure increase an oblique shock, although 
requiring a higher local Mach Number, has less total 
pressure loss than a normal shock. Consequently, 
unless the Mach Number is very nearly unity before 
the normal shock occurs, the total drag, including the 
effect of the relatively thicker boundary layer, will 
probably be greater with the combination of isentropic 
compression followed by a normal shock than it would 
be with the oblique shock alone. If the boundary 
layer cannot separate at all, the flow compression can 
be accomplished only by means of a large normal shock 
which would probably produce the greatest drag on 
the airfoil. However, the resultant drag, including the 
effect of the boundary layer separation, would probably 
be approximately of the same magnitude for all three 
cases and the type of flow would depend upon the 
boundary layer Reynolds Number. 

Since bs = C3; at Mach Numbers of 1.34 and 1.82, 
the oblique shock flow may be considered approxi- 
mately isentropic to the third order of flow deflection 
for Mach Numbers in this range. Fig. 12 for a Mach 
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Number of 1.5 shows that even near 6,,,, the pressure 
coefficient produced by an oblique shock is only 8 
*per cent less than that produced by isentropic com- 
pression. It is also seen that the fourth order approxi- 
mation from Eq. (32) is only 2 per cent less than the 
pressure coefficient produced by an oblique shock flow 
deflection of almost 10°. Fig. 13 for a flow deflection 
of 5° shows that the oblique shock and isentropic 
compression are in close agreement, and for Mach 
Numbers greater than 1.3 all the approximations of 
Eq. (32) higher than the first order give results prac- 
tically identical to that of the oblique shock. 

Fig. 14 shows that even for a +6,° value of 15° the 
isentropic compression will give a good numerical 
approximation for oblique shock flow for M; < 10. 
The power series approximation from Eq. (32) also 
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is good for Mach Numbers sufficiently greater than 
for Mach Numbers above 10 the 
third and higher orders of the power series representa- 
tion diverge rapidly from the oblique shock values. 
This occurs because the third and higher order odd 
terms of Eq. (32) have no limit as M, increases, as may 
be seen from Eq. (32) where for extremely large M, 


2 y+ (y + 1)? 
M ( 2 16 


37° — — 77 — 
6,4 


This shows that as long as 9, is finite, regardless of how 
small it may be, the third order term will increase 
without a limit as the Mach Number increases. There- 
fore, Eq. (32), Eq. (12), or similar expansion in terms 
of 6, cannot be used to investigate the flow at extremely 
large Mach Numbers. 


APPROXIMATE RELATIONS FOR EXTREMELY LARGE 
Macu NuMBERS 


Simple relations for the limiting values for an oblique 
shock as M, increases without limit may be obtained 


OBLIQUE SHOCK 


ISENTROPIC 
COMPRESSION 


15° 


+ 


PRESSURE COEFFICIENT, P 
T 


1 Do 


i 1 L x 
1.0 20 30 4050 10 20 30 4050 
MACH NUMBER IN FRONT OF SHOCK. mM 


Fie. 14. 


100 


from 


Ho 
that t 
tually 
flow 
heat 
intern 
formi1 
comp! 
rise W 
condu 
large 
an iso 
sumin; 
Eqs. | 
appro» 
with t 


Howev 
assume 
In this 
extrem 
large J 


or 


and fro 


It is see 
and (35 
P, nam 


as M, 


5 
\| 
| 
\ 
\ 
\ 
\ 
\ 
a 
| | | | 
; 9 
\ 
2 
SIN 8. 
+ 
. 
rf 


GE 


from Eqs. (20), (21), (24), and (27). 


Ps, (Qy/y + 1) sin® 
1 

M2? — (y — 1)/2y sin? — 
p2 tan 6,, | +1 (33) 


sin? 6, 
V2 


However, as M, increases without limit it is obvious 
that the physical limitations of the theory are even- 
tually exceeded. The so-called “‘exact’’ oblique shock 
flow equations were developed on the premise of no 
heat transfer (with the exception. of any occurring 
internally within the infinitesimal discontinuity sheet 
forming the shock wave). For large /, the resulting 
compression produces an extremely large temperature 
rise which would invalidate the supposition of no heat 
conduction or radiation. A better approximation for 
large M, would then be obtained by either assuming 
an isothermal process (equivalent to y = 1) or as- 
suming that 6, — @,. Either condition applied to 
Eqs. (20), (21), (24), and (27) yields the following 
approximation for large M, and finite @, which agree 
with those given by Epstein’ 

0.796, y = 1, Mo? 
Pt _, cin? — sin? 0, (Mi?) 
— Pr 
/opi Vi? 


However, for infinitesimal @, the process may still be 
assumed an irreversible adiabatic with no heat transfer. 
In this case, Eq. (28) or Fig. 6 shows that 6, becomes 
extremely small for small @,; therefore, Eq. (21) for 
large 4, becomes 


(34) 


— 2 sin? 0, — 2 sin? 6, 


A, — Pl 


or 
sin ~ — [(y + 1)/2]8, 


and from Eq. 


bo/bi + 1)/2]6,2Mi? 
Mz? — 2/y(y — 
(35) 
[(p2 — pr)/?/2p1 (vy + 


It is seen that for 0, + 0 and for y = 1 Eggs. (33), (34), 
and (35) give the same limit for the pressure coefficient 
P, namely: 


P = [(b2 — pi)/*/21 20,? 


as M, becomes extremely large. 
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In Fig. 14 the approximate relation for large M, is 
taken as 


P = (y + 1) sin? 6, = (12/5) sin? 6, (36) 


It is seen that this value closely approximates the 
“exact” oblique shock flow relation for M, > 15. 

It is interesting to note that for extremely high Mach 
Numbers Eq. (34) shows that the entire area affected 
by the oblique shock becomes only a thin film near 
the surface (6, — @,). Consequently, all the front 
surfaces of the body are exposed directly to M/; and 
Vi, and Eq. (34) may therefore be applied to all the 
front surfaces facing the stream. The pressure may 
be taken as zero for all the back faces and Eq. (34) then 
provides a complete pressure distribution for any two- 
or even three-dimensional body at extremely high 
Mach Numbers.’ 


PRESSURE DISTRIBUTION ON Two-DIMENSIONAL 
BopDIEs 


In supersonic flow the study of a two-dimensional 
body is even more important than in subsonic flow 
since the flow over any rectangular wing may be con- 
sidered everywhere two dimensional except at the 
wing tips. This is obvious because of the existence of 
the real characteristics in supersonic flow and the 
fact that nothing occurring downstream of a charac- 
teristic can directly affect the flow upstream from it. 
If the wing tip plan form is cut off along a Mach line 
then even the flow over the wing tip is two dimensional.® 
In addition, the flows over the upper and lower surfaces 
are now independent of each other because the condi- 
tions at the trailing edge can no longer directly affect 
the supersonic flow over the body. 

It is again worth while to emphasize that only the 
supersonic isentropic flow defined by Eqs. (3) and (5) 
is everywhere irrotational with a continuous velocity 
potential. The shock flow-defined by Eqs. (18) and 
(20) must always have a discontinuous drop in velocity 
potential at the plane discontinuity sheet defining 
the shock wave. The flow between shock waves can 
be irrotational with a continuous velocity potential 
only if the shock wave extends to infinity. However, 
for all practical purposes the shock wave for a finite 
body generally extends far enough from the body so 
that a continuous velocity potential may be assumed 
to exist in between the shock waves. This does show, 
however, that an exact rigorous solution is not yet 
available even when viscosity effects are neglected; 
consequently, the simplest approximate solution will 
probably be satisfactory within certain limits. The 
final criterion should, of course, be predicated on 
experimental verification. 

The following will attempt to provide a guide to the 
application of the solutions presented in this paper. 


(1) Osmaz.” > |0,°|< 2°, 1.1< Mi < 15 
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The Ackeret or linear approximation 


2 
P Cy WATE rs 1 6; 
will be satisfactory. This approximation is inde- 
pendent of the nature of the gas or the thermodynamic 
process and may be directly integrated to yield ex- 
tremely simple expressions for the forces on an airfoil.® 


(2) Osmar.° > 5°, 1.3< Mi< 15 
The Busemann or parabolic approximation 
P= CA, + C26,” = b0, + = 


2 + 1)/2] Mit — + 2 
6,? 
— 1 { (M,? — 1)? \ 


will be satisfactory. This second order approximation 
now depends upon the nature of the gas but is still 
independent of any entropy increase (irreversibility). 
This expression has also been formally integrated over 
the airfoil surface to show the dependence of the forces 
on the surface slopes.* ® 


(3) |@smar.” | > 6,” <= +15°, M, < 10 


In this case the third approximation given by the 
first three terms of Eqs. (12) and (32) may be used. 
This approximation now contains the effect due to the 
irreversible entropy increase caused by an oblique 
shock and is tedious to integrate formally. This has 
been done by Busemann* ® for a circular are profile. 
It may be seen that the value given by Busemann for 
D = (C; — bs) is incorrect as shown in Fig. 11 and the 
value of D should be taken from Eqs. (12) and (82). 
In addition, it is shown in Fig. 14 that the third and 
higher order approximations given by Eqs. (12) and 
(32) are invalid for M, > 10, or M; < 1.3. 

However, it has been shown that the flow approxi- 
mation given by Table 1 or Fig. 1 may be used for 
numerical evaluation of all the compression and ex- 
pansion flow over the body as long as Mi < 10, 0,° < 
15°. As previously pointed out, the use of Fig. 1 
automatically satisfies the conditions required for the 
existence of the solution. 


(4) |Osmaz.° | a 6,° > +15° 


Fig. 7 should be used to determine the pressure 
increase through all shock compressions where 0,° > 
15°, while Fig. 1 may be used for the regions of flow 
expansion and also for flow compressions where 0,° < 
15°. In this case the existence of the assumed shock 
wave and of supersonic flow throughout must always 
be checked from Table 2 and Fig. 9 although the use 
of Fig. 1 or Fig. 7 would automatically eliminate the 
need for Table 2. 

For cases (1), (2), and even (3) the flow is approxi- 
mately isentropic, therefore the pressure changes are 
reversible. Consequently, the pressure will always 
be the same when the surface slope (0,) is the same. 


For example, whenever 0, becomes zero, 2 = 1. 
However, for case (4) the process is now irreversible, 
and the pressure behind an oblique shock will be de- 
termined by Fig. 7 and cannot return to the corre- 
sponding value. For example, when @, again becomes 
zero behind a shock, p2 > fi. 


(5) M, > 15 


For large M, use Eq. (36) for flow compression and 
Fig. 1 for flow expansion. Then with extremely large 
My, use Eq. (34) for flow compression and for flow ex- 
pansion obtain a solution similar to Fig. 1 by utilizing 
the isothermal supersonic flow around a sharp corner.’ 
If the Mach Number is extremely large, then Eq. (34) 
may be used directly to obtain the entire pressure 
distribution for any two- or even three-dimensional 
body.’ 

Case (3) is the most important one for practical 
supersonic airfoil design. For the design of airfoils 
with sharp leading edges the use of Table 1 or Fig. 1 
may be greatly simplified if the variation of p/H, 
with @, is directly shown as in Fig. 15 and the pressure 
coefficient is taken as p/H rather than P. Then the 
movable @, scale may simply have its zero value shifted 
to the value corresponding to the free stream Mach 
Number (4) and the pressure coefficient (p/H;) is 
then read directly from Fig. 15 for the given surface 
slope. The inverse problem of determining the profile 
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that will produce any given physically possible pressure, 
distribution is also determined directly from Fig. 15. 
The scale for 6, in Fig. 15 either may be made movable 
or the difference in 0, required to shift the scale to a 
new M, may be algebraically added to the given surface 
slope. The position of the 6, scale is shown in Fig. 15 
for an M, of 1.5 and 1.7. 

For supersonic pressure distributions it is probably 
better to use the pressure coefficient p///; rather than 
P. The pressure summation to determine the resultant 
forces can be easily accomplished; for example, 


Mh = {1 + — = 


1 1 /2 M,? | 


therefore, 


(4) d (*) ins 
A, c 


Also from Eq. (11) 
P= 


11 + — — 1) 


The main advantage of the pressure coefficient P 
over p/H, is that a zero value of P shows the location 
of free stream values on the airfoil surface. However, 
in the type of supersonic flow assumed in Fig. 1 or 
Fig. 15 (0,° < 15°) the free stream conditions are 
always found whenever the airfoil surface has a 0, 
value of zero and are therefore already known. In 
addition, the local Mach Number over the airfoil 
surface is directly related to p/H by Eq. (9) or Fig. 1. 

Fig. 16 shows the pressure distribution obtained 
from Fig. 15 for the Stanton*® * Biconvex Circular Arc 
profile for a free stream Mach Number of 1.7. The 
p/H, value for the oblique shock at the sharp nose 
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/ 


(x/c¢ = 0) is obtained from Fig. 7 and also shown in 
Fig. 16. The p/H values in Fig. 15 are in good 
agreement with the '/.P values given by Hooker® 
except for a = 0°. Then the value at the nose (x/c = 
0, 8,° = 16°) as given by Hooker is much greater than 
that produced by an oblique shock. However, Fig. 15 
shows that even for this 0, value of 16° the p//; given 
by Fig. 14 is practically the same as for an oblique 
shock. Consequently, the values given by Fig. 14 
should have even better agreement with the experi- 
mental values than those given by Hooker. When 
Fig. 1 or Fig. 14 is used for 0,° > 15° the result should 
be checked by Fig. 7 for the p2/p: value and by Fig. 9 
to be certain that M; > 1. After the pressure forces 
on the airfoil itself have been obtained, the velocity 
and pressure in the wake may be approximated by 
utilizing Fig. 10 to obtain the total pressure loss through 
the oblique shocks. However, it must be noted that 
the boundary layer viscosity effects have been entirely 
neglected. 


If the maximum flow deflection (@mar.) for a given free stream Mach Number (M; = canst.) is desired, then 


Eq. (28) may be differentiated with respect to 6,, considering M, a constant and d0,/d0,, = 0. 


sin 20, — sin 


COS COS (I — 85) + Sin Sin (A. — | (37) 


cos? (6. — 


However, this is identical to the relation that would be obtained if the minimum free stream Mach Number 
(Mimin.) for the existence of an oblique shock for a given positive surface deflection (0, = const.) were required. 


Since for Mini». 1/Mi? is a maximum, Eq. (28) yields 


= (0/20.){sin? 0, — [(y + 1)/2][sin 6, sin 0,/cos (0. — = 0 (38) 
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Then both Eqs. (37) and (38) become 


y¥+1 sin 26, 
4 cos? (@, — 


[sin 26, — =0 (39) 


Eqs. (28) and (39) are then solved in the form 


sin 26. cos? (. Osmaz.) = 5 sin 20smar. 
6 sin 0,‘Sin Osmar. (40) 


5 cos (0. — Osmaz.) 


= sin? 6, — 
M 
and the values are given in Table 2 and Fig. 17. It is 
obvious physically that the solution of M,,,;,. aNd Osmaz. 
are equivalent when it is noted that Fig. 6 may be 
obtained directly from Eq. (28). 
An approximation to this exact solution was given 
by Busemann® as 
(4/3) (Minin. 
= 
3/2 
0.321 (Minin.* — 1) 


lmin- 


This approximates Table 2 for 0, values less than 20°. 
However, a more useful approximation may be obtained 
by observing in Table 2 that M; is only slightly less 
than unity for small @,. Consequently, the values of 
Table 1 approximate those of Table 2 for small @,, 
since Table 1 indicates the isentropic flow compression 
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Appendix II 


PowER SERIES APPROXIMATION FOR AN OBLIQUE SHOCK 


that would be required to reduce the Mach Number 
to exactly unity. Fig. 17 shows that the @,,,,.. ob- 
tained from the isentropic flow approximation (Table 1) 
agrees fairly well with the exact solution (Table 2) up 
tO a Dsnar. Of 20°. For Osmaz. less than 15° the isentropic 
compression is a desirable conservative approximation. 
Fig. 9 immediately shows that the isentropic approxi- 
mation is actually a better criterion for O@sn¢7. OF Mimin. 
in supersonic flow computations*with 6, < 15°. In 
Fig. 9 it is seen that for M; near M,,,,,. the flow behind 
the shock is still subsonic (Mz < 1); consequently, 
the oblique shock flow solution is no longer valid 
behind the shock. The use of Table 1 or Fig. 1 for 
Psmar. OF Mimsy. eliminates this difficulty by automati- 
cally keeping M2 greater than unity. 
The approximate relation 


Oy 1/4 + 
as given by Epstein’ may be obtained as an exact value 


for the limit of Eq. (40) as M, approaches infinity. 
When M, — ©, Eq. (40) becomes 


sin 20smar. Cos (A. — 9s) 
Sin Osmar. SiN By 
tan Osmar. = —COt 
therefore 
= + Omas./2 


However, as shown in Table 2, this relation is not even 
a rough approximation when M; < 4. As a matter 
of fact, Eq. (41) is not even physically justified for 
extremely large M, since it was previously pointed 
out that the “‘exact’’ oblique shock flow equations are 
invalid for large M, because the heat transfer accom- 
panying the extremely high temperatures produced by 
the shock is not considered in the theory. Actually, 
for large M, it would probably be best to use the 
relation Eq. (35) 


— [Cy + 1)/2]Osmac. 


and as a limit Eq. (34) yields 


lim 6, — 6; 


‘The pressure coefficient (P) may be obtained from Eq. (25) as 


P= (De Vi? = 2/[1 + (cot 6./tan 


Then from Eq. (27) 


sin? = {[(y + 1)/4]MvP + 1}/M? 


cot? = {(Mi?— 1) — [(y + 1)/4]M2P}/{[(y + + 1} 


(25’) 
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Therefore, assuming P is small and expanding in powers of P: 


cot = > a,P*— (42) 


n=1 


ai = — 1)'”, ag = —[(y + 1)/2°] — 1)'”7] 
(My? — 1)3/2’ 910 (M,? — 1)5/2 
_ (y + 1)4 (85Mi"* — 120M," + 144M," — 64M,") 
915 (M2 — 1)7/2 
+ (63M, — 280M, + 480M,"* — + 128M;j!*) 
gu on 1)9 /2 


a = 


Then, assuming @, is also small 


tan 0, = 0, + (0,3/3) + (20,5/15) + ..... 
and combining this with Eqs. (25’) and (42) gives 


P = 5,0," 
n=1 


by = 2/ay, by = —(b1/as)(1 + 
bs = —(1/a1)(b + 2asbibs + ashi? */s) 
by = —(1/ay)[*/sb1 + + + be”) + + 
bs = —(1/a1)[*/abe + by + + bobs) + + bo?) + + — (4/15)] 


bg = + 1/3b3 + + a2(2bibs + + bs”) + a3(6bib2b3 + + + + 602") + 
+ 


The coefficients 5; to bs could be computed for any M; directly from Eqs. (42) and (43). However, the coeffi- 
cients 6; to bs may be more conveniently written by combining Eqs. (42) and (43) to the form given in Eq. (32). 


It is also useful to write Eq. (43) in the form 


by b;3 


or from Eq. (32) 


(Mt 1)" + Mi + 1} p, 
2 4(M,? — 
Eq. (44) may be used to determine the angle of stream tube deflection if the local pressure coefficient P is known. 
Similarly, Eq. (42) may be transformed to an expression that would determine the local pressure coefficient if 
M, and 6, are measured, 


6; 


— 
P= (cot — VM? -l1)- (cot 0. — VM? — (cot My? | 
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INTRODUCTION 


Caw PROBLEM of determining bending moments in 
rotor blades is a fundamental one in the design of 
rotating wing aircraft. Several methods for solving the 
problem have been proposed, both empirical and 
theoretical. In general, the theoretical methods have 
been so involved that most blade designers have had to 
take recourse to empirical or semiempirical methods. 
Recently, Owen! has shown that, in several cases for 
which he obtained more or less precise solutions of the 
differential equation for blade bending, Cierva’s em- 
pirical formula'* gave excellent agreement with the 
“exact” solutions for the bending moment. Owen’s 
method is not, however, adaptable to the practical case 
of a blade of varying section. This has made it ex- 
tremely worth while to determine the assumptions and 
limitations inherent in Cierva’s formula, in order that 
it may be used with confidence in its range of applica- 
bility. 

It will be the purpose of this paper to review the 
fundamental equations that govern blade bending; to 
indicate several possible methods of solution, together 
with the simplifying assumptions necessary for each; 
and to discuss the analytical basis for various approxi- 
mate and empirical methods, particularly Cierva’s 
formula. Finally, an extension to hinged rotor blades 
of the Stuart-Myklestad tabular method of propeller 
blade stress analysis® '° will be made, together with 
some notes on numerical calculation. 

The question of aerodynamic loading on the blade 
will not be discussed here, since it has been treated in 
considerable detail in a number of other papers.® ° 


SYMBOLS 


EI = bending stiffness 
L = differential operator 


dr? dr? dr? 
M = moment on blade 
m = mass per unit length 
pb = axial loading (centrifugal force) per unit length 
Q, = total centrifugal force acting on the cross-sectional area 


atr 
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The Bending of Rotor Blades 
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radial distance from center of rotation to point on blade 

shear on blade 

transverse loading per unit length on blade 

transverse aerodynamic loading per unit length on blade 

inertia loading per unit length on blade 

rate of change of beam shear dS/dr 

blade bending deflection measured from line through 
hinge inclined at angle 8 

blade deflection measured from plane of rotation 

blade flapping angle 

blade azimuth angle measured from most downwind 

position 


ll 


ll 


THE FUNDAMENTAL EQUATIONS FOR BLADE BENDING 


The fundamental equations of blade bending may be 
derived by a consideration of Fig. 1. Here the blade 
resultant tranverse load per unit length is represented 
by ¢ and the axial centrifugal load is represented by p. 
In general, ¢ will be a function of r and y, while p will be 
a function of r only, for small displacements of the 
blade. The moment on a section through the blade at 
radius r is: 


M = — — S®pl(n — + (nn — dn 
(1) 
In accordance with usual beam theory, WM = EI(d*y/ 


dr*), or 


d*y R R 
= — r)dry — Pl(n 


Differentiating Eq. (2) twice with respect to 7, 
a? d*y f 


This same equation may also be derived by considera- 
tion of the rate of change of shear on a differential length 
of blade, as shown in Fig. 2. 
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The transverse loading ¢ is composed of two parts: 
one, t;, the aerodynamic surface forces, and the other, 


t,, the body forces, due to weight and inertia. Thus: 
(4) 
dt? dt? 


For blades of the type generally used in flapping rotors, 
the weight term, mg, is negligible in comparison with the 
aerodynamic and inertia loadings. The inertia term 
m(d*y/dt*) is in most cases small compared to mr X 
(d?8/dt*) unless the periodic aerodynamic forces are in 
resonance with one of the natural bending frequencies 
of the rotating blade. In this connection, it must be 
noted that the aerodynamic loading, 4, will also be 
altered somewhat by the deflection of the blade. At 
any section, the blade will have an additional transverse 
velocity dy/dt, which will change the local angle of at- 
tack and, hence, the aerodynamic loads. From the 
standpoint of resonance, this additional agrodynamic 
load will constitute a damping force. For most rotors 
conforming to current design practice, however, it ap- 
pears that y is small and resonance in the flapping plane 
does not occur, at least below blade stall. Measure- 
ments on the British C 30 autogyro? have confirmed 
this. 

Hufton, Nutt, Bigg, and Beavan? have included 
the effect of the d*y/dt? term in an analysis of the C 30 
rotor, but it does not appear to be important. In the 
analysis discussed in this paper, unless otherwise noted, 
the effect of blade flexibility on the aerodynamic and 
inertia loadings will be assumed to be neglected. In 
short, the transverse loading, 4, will be that calculated 
for a hinged rigid blade. 

The effects of torsional flexibility also appear to have 
been neglected in most blade bending analyses and will 
not be treated here. The influence of such torsional 
flexibility on blade loads, treated as a problem separate 
and distinct from the blade bending, may be found in 
papers by Wheatley* and Beavan and Lock.‘ 

In the solutions to the blade bending problem of 
Owen,! Hufton, et al.,? and Cierva,!* it is also assumed 
that the value of @ is that obtained for a rigid blade, 
which again merely amounts to assuming that y is 
small compared to Sr over most of the blade. In the 
Stuart-Myklestad® ® method, it is not necessary to 
make this assumption. This matter is, however, 
worthy of some further discussion, since it suggests a 
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means of simplifying the differential equation of the 
problem. 

The value of 8 is determined from the condition that, 
for a hinged blade, the moments about the hinge must 
be zero. Writing the moment at the root from Eq. (1) 
(r = 0), 


R R 
M= tr; dry — b (nB + n)dn = 0 (5) 


For a rigid blade y; = 0 or 


R R 
f pridr, = 0 (6) 


It may readily be shown that the omission of the 4; 
term in Eq. (5) is equivalent of the omission of the 
dy/dr term in Eq. (3). In the solutions of Owen, Huf- 
ton, et al., and Cierva, the value of 8 isdetermined on the 
basis of a rigid blade [Eq. (6)], but the dy/dr term is 
retained in Eq. (3). The dy/dr term is usually small 
compared to 8 (about 5 per cent maximum in a typical 
problem) so that it makes little difference in the solu- 
tion. When § is computed from Eq. (6), there appears 
to be little theoretical justification for either procedure. 
It will later be found convenient for certain analyses to 
neglect the dy/dr term. 

The basic equation to be solved in the theory of blade 
bending has been developed as 


ad? d*y 
(zt 0, + po pp (7) 
where Q, = JS;*pdr, is the total centrifugal force at 
any section. A mathematical solution of this equation 
for the general case where Q,, EJ, and p are arbitrary 
functions of 7 is extremely difficult; in fact, no such 
solution appears to have been attained even for the 
particular case of constant EJ and mass distribution. 
Hufton, et al.,? have pointed out the slowness of con- 
vergence of a series solution. Various approximate 
methods of analysis, most of them laborous, have been 
proposed. It is next proposed to review a number of 
these methods, together with their basic assumptions 
and limitations. 


VARIOUS METHODS OF SOLUTION 


Collocation 


Hufton, et al.,? have employed the method of colloca- 
tion in the solution of the problem. Basically, this 
method assumes that the solution is given by the 
finite series of n terms. 


y = Aigi(r) + Aada(r) + Anda(r) (8) 


where the ¢’s are arbitrary functions of r except. that 
they satisfy the boundary conditions and where the 
A’s are undetermined constants. Values of ¢,(r) satis- 
fying the boundary conditions may readily be formu- 
lated for most physical problems as algebraic poly- 
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5 nomials or trigonometric series. This expression for y ary conditions and ¢ is an infinitesimal constant. This, exc 
og 2 is then substituted back into the differential Eq. (7), of course, is nothing more than the mathematical formu- equ 
Rs leading to an algebraic equation with the mA’s as un- lation of the principle of virtual work. Using the mir 
= knowns. Then values of r corresponding to m arbi- Rayleigh-Ritz method for obtaining an approximate onl 
; trary points on the blade are substituted successively solution, y is assumed to be given by the finite series of of t 
a into the equation leading to equations in the m un- m terms so t 
‘4 knowns. It is thus assumed that the differential equa- c 
tion is satisfied at the boundaries and at the m arbi- = Aidilr) + Astalr) +... Anal?) of 
tos trarily chosen points on the blade. The degree of ap- where the ¢’s are arbitrary functions satisfying the obt 
a proximation is improved as more terms in the series boundary conditions and the A’s are constants to be 
, are chosen. Following this procedure, one is left with determined. Eq. (10) is then replaced by the condi- 
the unpleasant task of solving a considerable number of _ tion whe 
simultaneous equations. m=n tior 
In the solution by Hufton, et al., the terms in | EI (3) dr + boli 
dy/dr and d*y/d¢? are retained in the equation, but Bis 4» (0Am 0 2 0 2 \dr mu 
assumed to be the value obtained from the rigid blade R mat 
equations. dr — t'y ar =0 (12) tok 
its \ 
Type Solutions ‘ Upon substitution of Eq. (11) for y in Eq. (12),  simul- 
: Owen! has solved the problem for a blade of constant taneous algebraic equations in the #A’s result. a 
:s EI and mass distribution by an essentially inverse (Ggjerkin’s Method My 
method. Here a moment curve of some plausible sha 5. 
: is assumed, and the loading that ‘pbaithieth ‘nth Galerkin’s method’ is a generalization of the virtual —y 
5 moment is determined. Mathematically, this isa much work—Ray leigh-Ritz method for solving the differen- Den 
simpler problem, and Owen has determined related tial Am of mechanics. Galerkin , method is, 
loadings for a large class of useful moment curves, however, applicable to differential equations, whether 
which he has named “type solutions.” The procedure they arise from mechanics or not. Let us write Eq. (7) E] 
in solving a given problem is to find a type loading curve, in the form 
or a linear combination of such curves, which represents L(y) — t’ =0 (13) EI 


a satisfactory approximation to the given loading. where L is the differential operator 
The moments will then be given by the corresponding 2 Z @ d 
( ) dr? dr 


type moment curve or linear combination of such curves. pads 

The method is extremely useful for studying qualita- dr?\ dr? 
tive effects of changes in stiffness and blade weight for If an approximation for y other than‘the true solution 
blades of constant section, but it is not at all suited for is used, Eq. (13) will not be satisfied for all values of r. 
application to blades with arbitrary distribution of mass Instead, 


a and stiffness such as frequently arise in helicopter prac- 
tice. —t =e (14) 
where ¢ is the error of the approximation. The best Thu 
Energy Solutions : approximate solution is obviously the one that makes rect 
The method of virtual work may be applied to the the square of the error, integrated over the blade span, o_, 
problem. It can be shown by Euler’s equations of the ee ae 
Calculus of Variations that the exact solution of the 
: differential Eq. (7) is given by the function y, which a Liy) - | dr= ..* dr = minimum (15) yo 
annuls the first order variation of the integral expres- 
“a sion If the finite series of Eq. (11) is substituted for y, this Cier 
" condition may be expressed as 
EI (d*y\? RO. (dy\? Ci 
I= — dr + —{—) dr — R 
o o 2 \dr ra) rapic 
val L(y) - =0 (16) 
l tion 


This again leads to m simultaneous equations for the 
where t’ = ¢ — pB. This may be expressed as n unknown A’s. 


6I = 0 (10) Hohenemser’s Method 


when the variation is expressed by y’ = y + en(r), Hohenemser® has applied what is basically Galerkin’s 
where »(r) is an arbitrary function satisfying the bound- method to the solution of the blade bending problem, 
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except that he converted the differential equation to an 
equivalent integral equation before proceeding with the 
minimization of the error. In Hohenemser’s solution, 
only a single term is used for the unknown parameter 
of the problem, instead of a series similar to Eq. (11), 
so that, in general, close approximation is not to be ex- 
pected over much of the blade. For the case of a blade 
of constant EI and mass distribution, Hohenemser has 
obtained the approximate solution 


M = + 0.052(Q.R?/EI)] (17) 


where Mpg is the moment on a rigid blade. This solu- 
tion has been obtained for constant, linear, and para- 
bolic distributions of air load on the blade, without 
much change of the constant 0.052. It should be re- 
marked here, that Hohenemser does not assume 8 
to be the value for the rigid blade but rather determines 
its value for the flexible blade. 


Successive Approximations 


A method of successive approximations may be used 
to solve the integro-differential equation, Eq. (2), 
employing either analytical or numerical procedures. 
Denoting the successive values of y as, y, y™, y®, 

.. 9”, the following set of equations may be written 


dy R 
EI = t(r, — r)dn — — r)Bldn 


dy) R R 
EI = t(n— r)drn — ra —r) X 
B+ — y)]dr 


B+ — (18) 


Thus each succesive value of y may be obtained by di- 
rect integration. The process is continued until a 
satisfactory degree of agreement between y~! and 
y™ is obtained. Convergence will be rather slow be- 
cause the moment relieving effect of the centrifugal 
tension, which is neglected in the determination of 
y, plays so large a part in the final result. 


Cierva’s Method 


Cierva’s method'* was presented without proof as a 
rapid engineering method for obtaining rotor blade 
bending moments. The method is simple in applica- 
tion and consists of obtaining the bending moment for 
loads on a rigid but hinged blade, Mz, and the bending 
moment on a blade that is bent to the curvature it 
would have if only the centrifugal force resisted bend- 
ing, Vy Mp is readily obtained from Eq. (1) by set- 
ting y equal to zero. M, is obtained by the following 
procedure. EJ is set equal to zero in Eq. (7). 
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2+ pp (19) 
This is equivalent to 

(d/dr)[—Q,(dy/dr)] = t — pp 
or integrating, 

Q,(dy/dr) = (t — pB)dr (20) 


The integral on the left-hand side of Eq. (20) is the 
shear on a rigid blade as obtained from conventional 
beam theory. This may be denoted as S,. Then 


dy/dr = S,/Q, (21) 
The value of dy/dr is easily obtained, since both S, 


and Q, are basic quantities used in blade analysis. The 
value of M, is obtained as 
2. 
M, = = EI (22) 
dr? dr\Q, 


The differentiation of (S,/Q,) is usually done graphic- 
ally. Having the values of Mp and M, Cierva gives the 
following formula for the bending moment on the 
blade: 


M = + (23) 


The Stuart-Myklestad Method 


A tabular numerical method for propeller blade stress 
analysis has been presented by Stuart.? The method 
as presented originally required several repeated trials. 
Myklestad,”” however, has shown how the procedure 
may be modified so that the necessity for repeated 
trials may be eliminated. This method, although in- 
tended for application to propellers, is directly applic- 
able to blades having fixed root ends. It may be ex- 
tended to blades having hinged ends without much 
difficulty. 

In carrying out an analysis by this method, the 
blade is first divided into numbered stations at interval 
Ar apart. Midway between these stations are lettered 
stations as shown in Fig. 3. Then, from beam theory, 
to a degree of approximation which improves as Ar is 


made smaller: 
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2, — 2 = Ar(dz/dr), (25) 


Finally, the difference in moment between two stations 


1s 
Mz = M, + ArS, + (22 — 21)Q (26) 


The calculations may conveniently be set up in the 
tabular form given by Stuart. Dealing first with a fixed 
end blade, the moment ai all points on the blade may 
be determined if the root end moment is known. This, 
of course, is not the case. The tip moment, however, 
mus be zero, and any value of root moment other than 
the correct one will not give this result. In the Mykle- 
stad method, an arbitrary value of root moment, Mo", 
is assumed, and the calculations of Eqs. (24) to (26) 
are made with the air-load shears acting on the blade; 
the resulting tip moment is called Then asecond 
arbitrary moment, M,, is assumed at the root, and 
the calculation is repeated without any shears acting on 
the blade; the resulting tip moment is called M,®. 
The moment, M;, at any point on the fixed end blade 
is then given as 


My = M® — (M7) 


To take care of a blade hinged at the root, a unit 
value of dz/dr or coning angle is assumed at the root; 
no air-load shear is assumed to be acting; and a third 
set of moments on the blade, M“’, is obtained with 
the arbitrary root moment, M,“, and the arbitrary tip 
moment, M,. Since the tip moment must be zero 
in this case too, 


M® = M® ( M,* / M® 


where M“™ is the moment at any point because of the 
unit coning angle at the root. Now the moments on a 
hinged blade may be obtained from the condition that 
the value of root moment must be zero, which gives 


My = Mp — (My,/Mo)M“ 


In making the calculations by means of this method, 
it is necessary to use a considerable number of signifi- 
cant figures and to retain the same number of decimal 
places throughout the calculation. Five decimal places 
have been found to be satisfactory in a number of cases. 
In any particular solution, it appears to be desirable to 
run a check table on the final solution wherein the 
moments obtained are checked against blade curvature. 
Often it has been found that the values of moment check 
very well up to a certain station, and then the errors 
become increasingly large up to the tip where, of course, 
the check must again be satisfactory, since a zero value 
is assured by the method of calculation. When it is 
found that the moments do not check the curvature, it 
will usually be found satisfactory to repeat the calcu- 
lations using as the new assumed values of My"! and 
M,* the values obtained for My, and My. Also, the 
calculation may be repeated employing a larger number 
of decimal places. It has been found that these latter 
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procedures are more likely to be required with step 
tapered spars than with spars having no sudden changes 
of section. 


SOLUTION OF THE BENDING PROBLEM FOR A SPECIAL 
CASE 


For the special case where the ratio of centrifugal 
force to bending stiffness is constant, it is possible to 
obtain a simple solution for the bending equation, 
This case is of considerable interest for two reasons: 
First, it is a better approximation to actual conditions 
for a modern blade with step tapered spar than the 
common assumption of constant section; second, it 
provides a valuable insight into the theoretical basis for 
Cievra’s method, as will be shown later. 

If the dy/dr term is omitted from Eq. (7), (a step 
that is justified in the earlier part of this paper) and if 
t—p8 is replaced by w(r), there is obtained: 


a? d*y 

—| = 30 

From fundamental beam theory, 


M = EI(d*y/dr?) 
which when substituted in Eq. (30) leads to 


@M/dr*? — (Q,/EI,)M = w(r) (31) 
Let Q,/EI, be some constant, say k?; then, 
d?M/dr? — k?M = w(r) (32) 


The boundary conditions to be satisfied are 
M=0 atr=0 
M=0 atr=R 


The loading w(r) is arbitrary, but it is always zero 
at the root and at the tip, so that it may be represented 
by the following Fourier series: 


w(r) = A, sin (nxr/R) (33) 


where 


A, = — 34 
zh (w)r sin p r (34) 


in accordance with Fourier’s Theorem.'* A particular 
solution of Eq. (32) in operator form is 


R 


n=1 


The 
from 


Fro1 


or 


Usin 


App. 


This is easily solved to give 


+ k? R 


The complete solution to Eq. (32) is 
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(36) 


M=A cosh kr + B sinh kr + (r/R +B ™ 
sin (nmr/R) (37) 
By consideration of the boundary conditions, it is 
found that A = B = Oor 


nrr 


OR 


The bending moment on a rigid blade may be obtained 
from Eq. (32) by setting k = 0, which leads to: 


A . mmr 
Me = C, (39 
2. 


From the boundary conditions C; = C, = 0, and 


ol) 
= sin (mmr/R) (40) 


The bending moment on an entirely flexible blade, M,, 
as described by Cierva may be obtained from Egs. 
(30) and (22). 


M, = EI(d*y/dr?) = —EI(w(r)/Q,) 


(38) 


or 


M, = —w(r)/k? (41) 
Using the Fourier series for w(r), 
M, = — (A,/k*) sin (nr/R) (42) 
n=1 
Applying Cierva’s formula, 
A, . An . nur 
sin — sin — 
(nx/R)? R k? R 
M = (43) 


n=1 


This formula will be exactly equivalent to the “exact” 
solution Eq. (38) if the loading is given by a single har- 
monic of the Fourier series, or if Cierva’s method is 
applied separately to the moments resulting from each 
harmonic loading term, and the results added. In this 
connection it should be noted that for conditions that 
lead to maximum blade moments (aerodynamic and 
centrifugal components) the loading curve does not 
differ too greatly from a sinusoidal shape over the first 
half of the blade. 

Eq. (38) provides a fairly simple solution to the bend- 
ing problem for the blade of constant Q,/EI, or for 
blades that approximate this. The Fourier loading 
coefficients A, are not difficult to obtain using Eq. (34), 
if only a few terms are used, and, if greater accuracy is 
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desired, more terms may readily be obtained by some 
numerical scheme such as Runge’s method. The main 
value of this solution, however, lies in the insight into 
Cievra’s method which it provides. 


Analytical Basis for Cierva’s Method 


In order to provide a sound basis for the general 
comprehension and application of Cierva’s method, the - 
establishment of a theoretical basis for it seems es- 
sential. The results obtained for the blade of constant 
Q,/EI, may, in fact, be generalized so that no limita- 
tion need be placed on either Q, or EJ, except that they 
remain bounded and piecewise smooth functions, which 
is more a mathematical, than a practical, restriction. 

Suppose the variations in EJ, and Q, are given as 


Then 
= = (43a) 
Eq. (31) for blade bending may then be written 
(d?M/dr*) — yuq(r)M = w(r) (44) 
This is still subject to the boundary conditions 
M(0) =0 
M(R) = 0 
Now consider the equation 
(d?M/dr®) + yq(r) M = 0 (45) 


This equation for the specified boundary conditions can 
only have a solution for certain values of M, say 
Aa, A, AS is demonstrated in Courant-Hilbert.'! 
These values of u are called characteristic or “‘eigen”’ 
values. To each characteristic value \, corresponds a 
characteristic function ¢, which satisfies the equation. 

In the special case previously solved, where g is a 
constant, it may readily be seen that A, = (mr/R)? 
and 


= sin 


It can be demonstrated" that an equation of type 
similar to Eq. (45) will have an infinite denumerable 
number of \’s. Also, the system of functions ¢, will be 


orthogonal—i.e., 
= 0 (m 


The characteristics of the differential Eq. (45) dis- 
cussed above may be employed as indicated by von 
Karman (see reference 12) to obtain a solution of the 
equation. Returning to Eq. (44), w(r) is expanded asa 
series in the functions ¢, 


w(r)/q(r) = 


(46) 


(47) 


} 
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The coefficients A, may be determined by a method 
analogous to that used in determining the coefficients of 
a Fourier series, and as in the case of Fourier constants, 
the determination is made possible by the orthogonality 
condition, Eq. (46), plus the normalizing condition, 

= 1 (48) 


-Eq. (48) merely has the effect of fixing the amplitude 
of the @ function, so that the magnitude of the A,’s 
is set. Multiplying both sides of Eq. (47) by ¢2(7)q(r), 
the formula for A, is obtained as 


= (49) 


Thesolution to Eq. (44) may also be expanded as a 
series 


anda (50) 


Substituting Eqs. (47) and (50) into Eq. (44), 


n=1 n=1 


a(r) > (51) 


n=1 
Now, since each function ¢, satisfies Eq. (45) for u = Xn, 


: + = 


or 
= (52) 
Putting Eq. (52) into Eq. (51), 


which is satisfied if 
a, = + (53) 
Thus, 
M = ¥ + (54) 


A result entirely analogous to Eq. (38) is obtained for 


the special case. 
The rigid blade bending moment is readily obtained 


as 
Mp = — (55) 
The bending moment for the perfectly flexible blade is 
M, = (56) 


Now according to Cierva’s formula 


+ (An/u)]ba(0) 


Again, it may be seen that Cierva’s method agrees with 
the ‘‘exact’’ solution only if the loading can be described 
by a single characteristic function ¢, or if the loading is 
expressed as a series of such functions and Cierva’s 
method applied to each term separately. While it is 
possible to determine the ¢,(r) functions by a number of 
analytical and numerical methods, it is not suggested 
that this constitutes a practical method of analysis. 
The main reason for consideration of this method here 
is that it provides a mathematical foundation for 
Cierva’s method and indicates its inherent limitations. 

It may be concluded from our study’ that the accur- 
acy of Cierva’s method will depend upon the nature of 
the loading and upon the variation of the g(r) function 
for the blade. For the case of g(r) constant or constant 
ratio of centrifugal force to bending stiffness along the 
blade, which is of some practical value, the ¢ functions 
are of the form sin mzr/R. Where the value of g(r) 
does not vary greatly from a constant, ¢, should not 
vary greatly from the sinusoidal form. Thus, for a 
number of practical cases, Cierva’s method will become 
more or less accurate as the loading curve deviates less 
or more from a sinusoidal form. Further, Cierva 
moments determined over parts of the blade where 
sinusoidal restrictions on loading are complied with are 
in better agreement with the true bending moment than 
Cierva moments over other parts of the blade. 

What is actually done in applying Cierva’s method 
as indicated in Eq. (57) is to apply to the entire series 
an average correction of the type indicated in Eq. (54) 
for each term in the series of ‘‘eigen’”’ functions. If one 
term predominates, the Cierva method will automatic- 
ally give a correction nearly corresponding to this term. 
It is not necessary, of course, to know which term does 
predominate, or even what the ‘‘eigen’”’ functions are, 
in order to apply the Cierva method. If no one “eigen” 
function predominates in the loading, the error in ap- 
plying Cierva’s method will be somewhat greater, but 
the procedure will always lead to a mean correction. 
For this reason, it is much more satisfactory for ap- 
plication to arbitrary loadings than any single factor 
correction formula such as Hohenemser’s. The single 
factor correction formula generally corresonds to a 
single characteristic (or ‘‘eigen’’) function; if it is ap- 
plied to a loading in which other “‘eigen’’ functions pre- 
dominate, errors of larger magnitude will result. This 
may be readily realized by looking back at Eq. (38) for 
the blade of constant Q,/EI. If a single term correc- 
tion formula were derived on the basis of a first harmonic 
loading (n = 1) and then applied to an arbitrary load- 
ing in which the third harmonic (m = 3) predominated, 
the tremendous error is apparent. 
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constant section. 


COMPARISON OF VARIOUS METHODS 


Numerical computations were carried out by several 
methods in order to compare the accuracy of the 
method and the amount of labor involved. First, the 
blade of constant section chosen by Owen! for his illus- 
trative problem was investigated by comparing Owen’s 
solution with those obtained by the Cierva, Hohenem- 
ser, and Stuart-Myklestad methods. The results are 
shown in Fig 4 

Owen’s “exact” solution and the Stuart-Myklestad 
solution are found to be in good agreement. Owen’s 
results show slightly higher moments over most of the 
blade, particularly in the region of maximum moment, 
than those obtained by the Stuart-Myklestad method. 
This is due to the fact that Owen neglects the change in 
8 from the rigid blade value, whereas the Stuart-Mykle- 
stad method accounts for the change. 

The Cierva method compares favorably with the more 
precise methods over the first 65 per cent of the blade. 
It is considerably in error on the conservative side on 
the outer 35 per cent of the blade. Near the point of 
maximum moment, the Cierva method differs from the 
Stuart-Myklestad value by about 4 per cent on the 
unconservative side 

Hohenemser’s method shows the poorest agreement 
with the Stuart-Myklestad curve over most of the 
blade length and, in fact, does not indicate the correct 
point of maximum moment, which the other methods 
alldo. This is, of course, a result of the fact that the 
curve obtained by the Hohenemser method must be of 
the same shape as the rigid blade moment Surve. It is 


of interest to note that the Hohenemser curve differs Fic. 5. Blade bending moments for typical step-tapered spar. 
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from the more precise solutions first on the conservative 
side, then on the unconservative side, in such a way as 
to make the root mean square error small, which is in 
fact the basis of its derivation. 

In evaluating the merits of the various methods from 
a design standpoint, it should be remembered that with 
a flapping hinge, the design of the spar of constant 
section will be based on the maximum centrifugal force 
at the root or the combination of centrifugal force and 
bending moment at some point near the point of maxi- 
mum moment. Thus, the failure of Cierva’s method 
to predict correctly the moments on the outer part of 
the blade will be of little importance. It may then be 
concluded that, for blades of constant section and blade 
properties, Cierva’s method is a completely satisfactory 
engineering design method. Where greater accuracy is 
desired, the revised Stuart-Myklestad seems to offer 
the quickest “exact” solution and is more suited to 
manipulation by engineering personnel than the colloca- 
tion, energy methods, type solutions, or successive ap- 
proximations. 

For the case of typical step tapered spar, a compari- 
son of Cierva’s method with the Stuart-Myklestad 
solution for hovering is shown in Fig. 5. It may be 
seen that over most of the inner part of the blade the 
Cierva solution underestimates the moment by a 
maximum of 16 per cent in one spot. On the outer 
part of the blade, the Cierva moments are as much as 
250 per cent too high at r/R = 0.65 and still diverging. 
This is more serious in a step tapered spar than in a 
spar of constant section, since the outboard sections are 
designed largely by the value of the moment there. If 
weight is not too important, the spar could be designed 
by Cierva’s method allowing an appreciable margin over 
the inboard 50 per cent of blade and using the calcu- 
lated values over the outer 50 per cent. This is also a 
good basis for preliminary design. For an “exact” 
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solution, to obtain the minimum section, the re- 
vised Stuart-Myklestad method appears to be re- 


It should be kept in mind that the air-load data must 
be fairly reliable to justify accurate methods of blade 
design. At present, this is probably the greatest weak- 
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ness of blade design. Much work is being done on the 
subject, however, and it is hoped that improved air- 
load data will soon be available. 
engineer can only estimate the reliability of such air- 
load data from deflection and strain measurements on, 
and from observed performance of, past designs. 
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INTRODUCTION 


: wy AIRPLANE DESIGNER is faced with the selection 
of an airfoil section early in the design phases of 
each new airplane. The airfoil section is one of the air- 
plane design items which is seldom changed once an air- 
plane is built. This is necessary because so much of a 
modern airplane is intimately associated with the wing 
structure and changing the airfoil is usually tantamount 
to designing and building a new airplane. Therefore it is 
important that a reasonably good choice of airfoil is 
made in the shortest possible time at the start of each 
design. 

The designer finds that the wing on the average air- 
plane is probably the primary factor in defining the air- 
plane general arrangement. First of all the span, area, 
dihedral, sweep, and taper are important. Then the air- 
foil section itself is important. The section usually 
varies from root to tip and this variation is also im- 
portant. Some of the more significant airfoil factors 


(1) Drag. The drag of a wing is principally con- 
trolled by the airfoil section and its surface condition. 

(2) Lifts. The primary purpose of a wing is to 
produce lift. The lift characteristics usually control the 
stability. 

(3) Moments. The moment of the wing is struc- 
turally important and affects the stability. 

(4) Stalling. The airfoil section usually defines the 
stalling characteristics of an airplane. 

(5) Compressibility. The airfoil section is the 
primary factor determining the effect of Mach Number 
on airplane characteristics. 

(6) Flaps. The ability of an airfoil to be highly 
flapped is frequently important. 

(7) Ailerons. Aileron design is intimately associ- 
ated with airfoil design. 

(8) Structure. Airfoil proportions, particularly 
thickness, are important in defining structure. 

(9) Construction. Airfoil proportions greatly af- 
fect the simplicity and ease of construction. 

(10) Landing Gear. In many airplanes the wing is 
called upon to house the landing gear and this imposes 
space requirements upon the airfoil. 

(11) Wing Volume. The wing is usually called upon 
to house fuel and sometimes engines. The volume re- 
quirement will frequently design wing thickness. 
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(12) Deicing. The problem of wing deicing is 
directly related to leading edge shape. 

(13) Flutter. The flutter characteristics of -an air- 
plane are frequently controlled by airfoil proportions. 


TESTED AIRFOILS 


The airplane designer has almost a limitless number 
of tested airfoils to choose from. The National Ad- 
visory Committee for Aeronautics has been active since 
its beginning in presenting data on airfoils. Many of 
these airfoils were tested in systematic series and the 
results have been presented in reports such as N.A.C.A. 
T.R. Nos. 460 and 610. In view of the large number of 
airfoils tested it would seem that the field should be well 
covered and little room for improvement should re- 
main. Such a stand does not bear much investigation. 
There is still much to be learned about the behavior of 
airfoils. A great deal of theoretical work has been done 
to predict and explain the flow over airfoils. It has been 
found that theory and experiment agree very well in 
flows where separations do not exist. In addition, 
thick boundary layers cause divergence of experiment 
from theory. These theories, in general, are primarily 
for incompressible flows, and extensions to include the 
effects of compressible fluids are only satisfactory when 
no part of the flow approaches sonic conditions. How- 
ever, most airplanes are designed to operate without 
flow separation on the airfoil except when approaching 
the stall and are flown at speeds where compressibility 
effects are reasonably accounted for by theory. This 
means that, with usable airfoils, it is possible to predict 
by theory all of their characteristics except the stall 
and compressibility effects. A partially empirical pres- 
entation of the basic characteristics of common airfoils 
is given by Jacobs and Rhode in N.A.C.A. R.T. No. 631. 


ANALYSIS OF AIRFOIL DATA 


A typical set of airfoil data as measured in a two- 
dimensional wind tunnel is presented in Figs. 1 through 
7. These data will be analyzed to show how a few tests 
can be used to indicate the characteristics of many air- 
foils. The data that are sometimes available may in- 
clude: 

(1) Two-dimensional wind-tunnel data giving lift, 
drag, and moments at one or more Reynolds Numbers. 

(2) Pressure distributions for (1). 

(3) Two-dimensional wind-tunnel data giving lift, 
drag, and moments at various Mach Numbers. 

(4) Pressure distributions for (3). 
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Fic. 1. B.A.C. No. 53, airfoil and ordinates. 
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airfoil. 


(5) Effects of surface condition and turbulence on 
characteristics as measured in two-dimensional wind 
tunnels or in flight. 

(6) Three-dimensional test data on wings incorpo- 
rating the airfoil. 

The sample data presented in Figs. 1-6 do not in- 
clude all possible tests but are sufficiently complete to 
be useful in understanding some of the phenomena to 
be encountered. An analysis of the pressure data ex- 
plains many things. In many ways, pressure data are 
easier to use when the pressures at any point are plotted 
versus angle of attack or normal force.coefficient which 
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Cross- 


Fic. 3. Cz vs. Mach Numbers. B.A.C. No. 53 airfoil. 
plot of Fig. 2. 
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Fic. 4. Cp vs. Cy at various Mach Numbers. B.A.C. No. 8 
airfoil. Drags measured by wake survey. 


is usually even better. The data of Fig. 6 are replotted 
versus normal force coefficient in Fig. 7. It will be no- 
ticed that, for nearly all points, the variation of pressure 
coefficient with normal force coefficient is approximately 
linear. Theory, however, indicates that the velocity 
and not the pressure should be linear. The pressures 
can be converted to velocities by the use of Bernoulli’ 
equation, which indicates that the square root of the 
pressure coefficient will be the ratio of velocity to free 
stream velocity. The velocity ratio or square root d 
the pressure coefficient is plotted versus normal fore 
coefficient in Fig. 8. In this plot the variations at 
more nearly linear. Theory indicates, and experimett 
confirms, that the variation of velocity coefficient with 
lift is nearly identical at each per cent chord station for 
all airfoils. In general, differences occur only in tht 
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Cp vs. Mach Numbers. B.A.C. No. 53 airfoil. Cross- 
plot of Fig. 4. 


Fic. 5. 


first 1 or 2 per cent of airfoil chord. Since pressures are 
more useful experimentally, it is convenient to plot 
pressures on graph paper with a distorted scale so that 
a vertical scale is linear in velocity ratio rather than 
pressure. Examination of pressure data on a large 
number of airfoils indicates that appreciable variations 
in the pressures from those given by extrapolations us- 
ing standard slopes of velocity with lift occur only when 
separation begins. The pressures or velocities at zero 
lift are a function of the airfoil shape. If the airfoil were 
symmetrical, the pressures or velocities would be sym- 
metrical and, therefore, any differences from symmetry 
can be ascribed to the camber. From tests on a few air- 
foils it is possible to determine the velocities for many 
camber lines and thickness variations. For moderate 
variations in camber or thickness ordinates the veloci- 
ties can be assumed to be proportional to the ordinate 
when all ordinates are uniformly increased or de- 
creased. 
The characteristics of the airfoil without separation 
can be related to the physical shape by means of the 
pressure distributions. The slope of lift curve for infi- 
nite aspect ratio must be estimated from experimental 
data but is not susceptible to much variation. The mo- 
ment coefficient is constant about some aerodynamic 
center. The value of this moment coefficient will be in- 
dependent of aerodynamic center at zero lift and is 
therefore a function only of the camber. From the 
Pressure data the moment coefficient for each camber 
line can be computed. The aerodynamic center is de- 
fined by the variation of pressures with lift, and, since 
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Fic. 6. Pressure distribution over B.A.C. No. 53 airfoil. 
0.10, R.N. = 1.43 X 108, 
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Fic. 7. Surface pressure P/g vs. Cy for B.A.C. No. 53 airfoil. 
M = 0.10, R.N. = 1.43 X 106. 


this is nearly identical for all airfoils, there will be little 
variation of aerodynamic center between airfoils with 
the center usually lying just ahead of the quarter 
chord. The effect of camber on the angle of attack for 
zero lift is obtainable directly from the camber data. 
Drags cannot be computed so easily. When separation 
does not exist, the drags will be determined primarily 
by the relative amounts of laminar and turbulent flow. 
In some cases pressures will define the transition point, 
but in nearly all practical cases the transition will be de- 
termined by constructional details and by such items as 
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with full knowledge of the test Reynolds Number and 
possible effects of its variation. Much information on 
this subject has been given by the N.A.C.A. in reports 
such as T.R. No. 586. Some classes of airfoils are sus- 
ceptible to Reynolds Number effects, while others are 
not. In general, increasing the Reynolds Number in- 
creases the maximum lift and delays separation. The 
stall characteristics eventually deteriorate at high Rey- 
nolds Number. If all other things remain equal, it will 
usually be desirable to use an airfoil not greatly affected 
by Reynolds Number, since in this case wind-tunnel 
tests of a complete model can be expected to give reli- 
able stability and control data. Although drag sav- 
ings, indicated by model tests in regions where separa- 
tion exists, may not actually exist full scale, it will nearly 


9 ft., as are currently used for nearly all development 
testing. The converse is not always true. Model de- 
fects dc not always show up full scale, but to assume 
that they will not is dangerous. In most cases, the tuft 
appearance of separations and the stalls is not greatly 
affected by Reynolds Number. Usually, the difference 
between good and bad stall characteristics lies in the 
rate of progress and location of the separation. These 
can be easily evaluated in model tests and are impor- 
tant criteria to be examined. 


COMPRESSIBILITY 


Compressibility characteristics can be evaluated very 
well in most cases from two-dimensional tests. Fig. 9 
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Fic. 9. Comparison of two-dimensional and three-dimensional 
airfoil tests. B.A.C. No. 53 airfoil. 


compares two-dimensional with three-dimensional air- 
foil data for the B.A.C. No. 53 airfoil. The lift data 
are by forces in both cases. The drag data are by 
force in the three-dimensional case and by wake survey 
in the two-dimensional case. Excellent agreement has 
been found in several such cases. When the three-di- 
mensional test conditions are not favorable, better re- 
sults can be obtained two dimensionally. Such a case 
is shown in Fig. 10 where the three-dimensional model 
results were badly affected by blocking in at least one 
case. None of the data presented in this report has been 
corrected in any way for blocking. It is of great inter- 
est to be able to predict the compressibility character- 
istics. Using the relationship of von Karman for the 
critical Mach Number (Mach Number for attainment 
of supersonic flow), it is possible to compute Mach 
Numbers for which supersonic flow is just attained. 
Using the Karman-Tsien relationship to predict this 
critical Mach Number and the Glauert relationship to 
compute lift coefficient, it is possible to obtain estimates 
of the critical Mach Number and the lift coefficient at 
that Mach Number corresponding to a given low-speed 
pressure peak and lift coefficient. A convenient chart 
for making such calculations has been developed. For 
each low-speed pressure peak value there will be a 
critical Mach Number in accordance with the Karman- 
Tsien relationship. At this Mach Number the lift co- 
efficient is increased over the low-speed value by the 
Glauert relationship. 
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Fic. 10. Comparison of two-dimensional and three-dimensional 
airfoil tests. B.A.C. No. 100 airfoil. 


WING SECTION DEVELOPMENT 55 


In the application to an airplane, it is speeds and 
wing loadings that are of interest and not Mach Num- 
bers and lift coefficients. Mach Number can be easily 
converted to speed as accomplished at the right of Fig. 
8. By definition W/S = (p/2)V°?C,. This is readily 
converted to W/S = (p/2)a*M?C,. It can be shown 
that a*(p/2) = p(y/2). Then W/S = p(y/2)M°C, or 
W/S = (b/po) 1,485M?C,;. On Fig. 8 lines of constant 
M?C, have been drawn for conditions of lift coefficient 
and Mach Number at the start of supersonic flow. The 
lift coefficient scale at the bottom of the page is for low 
speeds as just explained. At the top of Fig. 8, scales are 
given to interpret °C, as wing loading at various alti- 
tudes. This chart is used by plotting upon it either the 
envelope of pressure data or the complete pressure data, 
as in Fig. 8 and explained previously. The graph has 
been made to have a linear velocity ratio vertically. 
From this plot a given wing loading and altitude may 
be chosen at the top and followed downward along a 
constant M°*C, line to the envelope of experimental 
data. Then, moving horizontally to the right, the 
critical Mach Number and speed at the altitude chosen 
can be read. In this manner the low-speed pressure data 
can be rapidly evaluated. It is interesting to note that 
in some cases the constant lift (17*C,) lines cross the 
airfoil data at two points. There is a high speed above 
which supersonic flow is obtained and a low speed below 
which it also occurs. The meaning of this is of great 
importance, for it can be seen that there is a wing load- 
ing above which supersonic flow will occur regardless of 
speed. 

This method permits the determination of speeds at 
which supersonic flow should occur. These critical 
speeds are not necessarily related to the speed at which 
the drag starts to increase or to the stability critical 
speed. Actually, the drag starts to rise at about the 
point where supersonic flow starts. This drag rise 
comes about not only because of the shock, where the 
flow attempts to become subsonic again, but also be- 
cause of separation in the flow. The pressure gradients 
are increased by Mach Number, and, if separation is 
imminent, this can cause a drag rise without the exist- 
ence of supersonic flows. If separations are not likely, 
there may be no drag rise at all up to the point where 
supersonic flow starts. The relationship between pre- 
dicted speeds and experimental values is shown on Fig. 
5 for the data of Fig. 8. 

The method for obtaining two-dimensional airfoil 
data at high Mach Numbers varies between various ex- 
perimenters. An unusual method has been used by 
the Boeing Aircraft Company in their wind tunnel at 
Seattle. Two walls, 6 in. apart, are mounted vertically 
in the test section giving a test section 6 in..in span by 
96 in. high. When 6-in. chord airfoils spanning the jet 
are tested in this tunnel, the blocking is small, per- 
mitting tests of 10 per cent thick airfoils to speeds 
above 0.92 Mach Number. Figs. 11 and 12 show the 
arrangement of this test section. The velocity be- 
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Fic. 11. High-speed two-dimensional end plates in test section 


of Boeing wind tunnel. 


Fic. 12. Top view of high-speed two-dimensional test section 
with tunnel top removed. 


tween the walls is uniform at all speeds and is about 50 
per cent higher between the walls than in the tunnel 
proper. Lift, drag, and moment are measured on the 
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force balance. The drag is also measured by wake sur. 
vey. Figs. 13 and 14 show the test section used to 
measure low-speed pressure distribution, flap effects, and 
stalling characteristics. Airfoils used in this test sec. 
tion have 36-in. span and usually are of 24-in. chord. 


APPLICATION TO THE AIRPLANE 


So far, only the determination of the aerodynamic 
characteristics of an airfoil have been covered. Meth- 
ods have been described by which the physical charac. 
teristics of an airfoil may be resolved into aerodynamic 
characteristics. However, the physical problems are of 
great importance in determining the airfoil. In most 
airplanes, the wing is called upon to be the backbone of 
the structure, to be the container of fuel and landing 
gear, and to be the hiding place for a high lift device, 
as well as to produce lift and lateral rolling control, 
The flexural rigidity of a wing is closely related to the 
ratio of span to root thickness. To keep equal rigidity 
when changing aspect ratio, it is necessary to thicken 
the root sections proportional to the aspect ratio in- 
crease. This leads to the use of high thickness ratios 
for the root sections of high aspect ratio wings. Most 
long-range airplanes have fuel volume requirements 
that are so severe that wing thickness is defined by the 
thickness necessary to provide the required internal vol- 
ume. If the airplane is not expected to fly at high Mach 
Numbers, the maximum wing thickness may still be 
limited by aerodynamic considerations. The maximum 
possible thickness ratio has not yet been determined. 
The Boeing XPBB-1 used a constant chord inner wing 
with a thickness of 23 per cent. Entirely normal flow 
existed over this wing even with the extremely large de- 
icers inflated. 


LANDING GEAR AND FLAPS 


The landing gear is frequently housed in the wing. 
This will often place severe limits upon the airfoil pro- 
file if the gear is to be held within the normal wing sur- 
face. Usually, a compromise is reached with portions 
of the landing gear protruding from the wing contour. 
The problem is more serious with wing flaps where the 
wing contour is usually chosen se as to contain all or 
nearly all the flap mechanism. Hollowness of either the 
upper or lower surface contours frequently works 
great hardships upon the flap mechanism. Such mech- 
anism problems will frequently be more important 
than aerodynamic considerations in defining the 
trailing-edge shape. 
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The aileron normally fits into the outer 40 per cent of 
the span of the trailing edge. This places additional 
demands upon the airfoil sections of this portion of the 
wing. Airfoil sections with large included angles at 
the trailing edge tend to have thick boundary layers 
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at the trailing edge with a strong tendency toward local 
separation. This will greatly affect the aileron ef- 
fectiveness and hinge moments. If the airfoil camber 
line has appreciable camber near the trailing edge, the 
aileron hinge moments are likely to be biased the same 
as they would by a full span tab with similar shape. 
This will cause the ailerons to float at other than neu- 
tral and will load the aileron control system. This is 
undesirable on military airplanes where battle damage 
may sever a control cable. It also causes added control 
system friction. When spring tab ailerons are used, it 
is important for the ailerons to float at neutral. The 
other aerodynamic effect of the airfoil on the aileron is 
the effect of bulging or hollow contours on hinge mo- 
ment. Optimum results are obtained with nearly 
flat surfaces. For all of these reasons, it appears desir- 
able to use an airfoil that has flat surfaces in the re- 
gion of the aileron with a fairly small included angle. 
The use of flat surfaces to form the trailing edge of 
airfoils has much to be said for it, although it places 
limitations upon the camber and thickness shapes that 
can be used, since both will automatically become 
straight in the region where flats are used. In general, 
such limitations are entirely acceptable. The effect 
of such a change in airfoils near the trailing edge cannot 
have any great effect upon drag because the boundary 
layer in this region tends to be flat even with curved 
airfoils. This has been borne out in two-dimensional 
airfoil tests conducted at Boeing. When the wing has 
a constant .thickness ratio, flat surfaces are useful in 
that the trailing edge becomes two flat planes. This 
permits simplified structural design and is particularly 
useful when flaps or ailerons are not made of a constant 
percent chord. In this case, it is possible to build con- 
stant chord flaps or ailerons with identical shapes along 
the span. 


FLUTTER 


The flutter and aileron reversal problem is intimately 
associated with the airfoil shape. The torsional rigidity 
that can be built into a wing is directly related to the 
cross section of the airfoil. Thin airfoils are fundamen- 
tally flexible and therefore undesirable from an aileron 
teversal standpoint. The flutter characteristics are also 
affected by the chordwise location of the structure in 
that this locates the mass center and the shear center to 
a great extent. This factor is not normally considered 
in the choice of an airfoil but may become important as 
speeds increase. 


DEICING 


Systematic studies of the effect of different airfoils 
upon deicing have not been made in deciding upon 
airfoils, but this may be important. Sharp leading edges 
should materially reduce the area to be deiced. There 
appears to be ample evidence to indicate that airfoil 
contour has much to do with the type of ice accumula- 
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important effects upon the airfoil choice. 
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Fic. 13. Test section of Boeing wind tunnel with end plates for 
low-speed tests of 2-ft. chord = 3-ft. span airfoils. 


Wake survey rake behind flapped airfoil in test section 
of Fig. 13. 


tions encountered. This problem should be studied in 
ice tunnels to determine the airfoil shape factors af- 
fecting icing. 


STRUCTURE 


The structural design of wings may frequently have 
In the sim- 
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plest case, this would result in the use of a rectangular 
wing with a flat lower surface (very common 10 years 
ago), but it is also likely to affect thickness distribution 
and surface curvature. Peculiar twists and nonuniform 
curvature will frequently be so difficult to form that the 
wing incorporating them may be so distorted when it 
comes out of the shop as to definitely indicate some sim- 
pier shape would have given better aerodynamic re- 
sults. This manufacturing problem must not be over- 
looked in choosing airfoils. 


STALLING 


Usually the first consideration in the choice of air- 
foils is the stalling problem. Stalling characteristics 
usually define the taper and twist of the wing for the de- 
signer. If a designer is prejudiced favorably or other- 
wise with respect to any particular airfoil, it is nearly 
always because of some experience with the stalling of 
an airplane using the airfoil. Although taper and twist 
are important, there can be little question that the air- 
foil is the most important single variable in the stalling 
problem. In multiengined airplanes, the propeller slip- 
streams will tend to prevent root stalling when any power 
is being used. The stall is thereby forced onto the tips, 
and the problem is not one of obtaining root stalling but 
one of having good stall characteristics with the stall first 
appearing at the tips. This is contrary to the popular 
conception wherein ‘“‘tip stalling’ is used to define bad 
stalling characteristics. In solving this problem, the 
airfoil shape is almost the only variable which can be 
varied. Airfoils with flat top lift curves give progres- 
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sive stalls with more than sufficient warning to 
pilot. The conventional approach to computing the 
stall by the use of spanwise lift curves has little if am 
significance when applied to wings with flat top li 
curves. The application of flaps to airfoils material 
changes their stall characteristics. This is likely to kk 
important in the choice of airfoils and flaps. In any 
case, the airplane designer will normally use only thos 
airfoils that for one reason or another he thinks wil 
give a good stalling airplane. Although the stalling 
problem is of probably first importance in the choice @ 
airfoils, it is the one concerning which the designer has 
probably the least factual information. 


CONCLUSIONS 


Each airplane and each designer will face different 
problems in the choice of airfoils. Some of the problems 
that are faced and some systematic methods for their 
study have been discussed. The properties of sections 
without supersonic flow or separation are readily pre. 
dicted. The critical Mach Number can be estimated 
reasonably well, but the boundaries and characteristics 
of separated flow are not yet susceptible of estimation. 
Although these separation problems define the stall, 
which is of primary importance, experimental methods 
for determining stalling characteristics are questionable. 
Physical, structural, and arrangement problems will 
often define airfoils and their thickness in particular. 
Aileron and flap problems can have strong effects tend- 
ing toward the use of airfoils with straight contours on 
their rear portions. 
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ABSTRACT 


It is now common practice to measure forces by a spring in the 
form of a cantilever beam whose deflection is measured by re- 
sistance wire strain gages. In this paper, formulas and charts 
are given to aid in designing the beams. Several electrical cir- 
cuits are analyzed. The results of some experiments on the ef- 
fect of working voltage and temperature are presented. % 


NOTATION 


F = applied force, lbs. 

W = beam width, in. 

Aa/a = unit strain in beam, in per in. 

E = Young’s modulus, lbs. per sq. in. 

h = half thickness of beam, in. 

% = distance from point of application of load, in. 
6 = deflection, in. 

M.A. = Mechanical advantage. of lever system 

R = resistance, ohms 

R,, Ro, Re gage resistances, ohms 


R;’, Ro’, R3’, Rs’ = lead resistances, ohms 

AR = change in resistance, ohms 

= strain sensitivity, AR/R = S(Aa/a) 
= output voltage, volts 

= input voltage, volts 

= meter resistance, ohms 

= output current, amp. 

= input current, amp. 

= shunt resistance, ohms 

temperature, °C. 

a’ = temperature coefficient, ohms per °C. 
= applied load, lbs. _ 

= length of wire in gage, in. 

= diameter of wire, in. 


ll 


INTRODUCTION 


es SEVERAL YEARS resistance wire strain gages have 
been extensively used for measuring stresses in 
structural members. The methods used have been ex- 
tensively covered in the literature. Recently, the air- 
craft industry and the National Advisory Committee 
for Aeronautics have been using resistance wire strain 
gages for measuring forces in both flight and model 
tests. These methods have now been sufficiently de- 
veloped to offer possibilities of use in secondary force 
measuring systems, such as wind-tunnel balance sys- 
tems and the like. In this faper it is proposed to pre- 
sent some of the principal considerations involved in 
such applications. 


Presented at the Radio and Instruments Session, Fourteenth 
Annual Meeting, I.A.S., New York, January 29-31, 1946. 
* Aeronautical Engineer. 
t Research Analyst. 
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APPLICATIONS 


Forces can be measured by springs whose deflection 
is measured by resistance wire strain gages. Two types 
of springs are commonly used—the cantilever beam and 
the ring. Of these, the cantilever beam offers more ef- 
fective use of its surface area and usually presents no 
installation difficulties. This discussion will be confined 
to cantilever beams. 

A constant strength beam of parabolic shape has the 
advantage that the average strain is equal to the maxi- 
mum strain. In practice, two points are chosen on the 
parabola and a straight line drawn between them. For 
design purpose, it is sufficient to use the formulas for a 
constant strength beam. Corrections have been cal- 
culated and expressed graphically but their use has not 
been found necessary. 

The handbook formula for stress in a constant strength 


beam is 


F/W = */3( Aa/a)E(h?/x) (1) 

The corresponding formula for deflection is 
8°(F/W) = (6/27)E(Aa/a)*x® (2) 
In 


These formulas are graphed in Figs. 1 and 2. 
practice the deflection and maximum load are usually 
specified and the length and width are set by practical 
considerations. The corresponding unit strain is de- 
termined from Fig. 1 and the beam thickness from Fig. 
2. The quantity Aa/a determines the sensitivity of 
the system. If this turns out to be too small, a suitable 
lever system will raise the sensitivity, for after a lever 
system is introduced, we have the new force and de- 
flection at the beam. Then from Eq. (2) 


F’ = F/M.A., 6’ = 6M.A. 


Aa\? 
a 


which shows that the sensitivity for a given load and 
deflection varies as the cube root of the mechanical ad- 
vantage. 

The most useful electrical circuits are based on the 
Wheatstone bridge. It is most convenient to have all 
four arms of the bridge active. This makes temperature 
compensation simpler and, in some circuits, gives more 
sensitivity for a given amount of energy dissipated on 
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the beam. Three circuits that have been used success- 
fully are illustrated in Figs. 3, 4, and 5. 

The circuit of Fig. 3 uses a potentiometer to measure 
the bridge unbalance. The sensitivity is given by 


e = s(Aae/a) 


When used with modern self-balancing potentiometers, 
this is usually the most satisfactory arrangement. The 
sensitivity of the circuit is independent of the tempera- 
ture. The zero drift with temperature, as will be seen 
later, is easily eliminated. If a potentiometer with a 
sensitivity of 0.01 millivolt is used with an applied volt- 
age of five volts and gages with a gage factor of two, this 
setup will measure a (Aa/a) = 1 X 10, which for 
steel corresponds to a stress of 28 Ibs. per sq.in. 
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The circuit of Fig. 4 uses a galvonometer or micro- 
ammeter to measure the bridge unbalance. The sensi- 
tivity is given by the formulas 


= {s(Aa/a)/[1 + (r/R)]}J 
e = {s(Aa/a)(r/R)/[1 + (r/R) }}e 
ei = {s(Aa/a)*(r/R)/[1 + (r/R)]*}el 


The last of these, the power sensitivity, is the most im- 
portant. The bridge and meter must be designed to 
go together. For maximum efficiency 7 and R should 
be equal. However, in this case, the circuit is usually 
overdamped, therefore any circuit of this type must 
be a compromise. With a commercial microammeter 
a sensitivity of Aa/a = 1 X 10-5 can easily be obtained. 
Much more sensitivity may be obtained with galvanom- 
eters. This circuit is used for measuring transient 
forces. At medium frequencies, the bridge output may 
be fed directly into a galvanometer. For high frequen- 
cies, amplifiers may be necessary because of the low 
sensitivity of high-frequency galvanometers. The 
sensitivity of the circuit depends on r/R and, hence, 
on the temperature; however, for the low temperature 
coefficient gages ordinarily used, this effect is negligible. 

The circuit of Fig. 5 uses a galvanometer or potentiom- 
eter only as an indicator and the shunting resistance is 
varied to give a zero reading on the indicator. The 
sensitivity is given by 

(Ap/p)(R/p) 
a@ [1 + (R/p)\(2 — {1/[1 + R/p]}) 


It has the advantage that the sensitivity is independent 
of the supply voltage and therefore a separate voltage 
reading is necessary. On the other hand, the calibra- 
tion curve will be nonlinear and the sensitivity again 
depends on the temperature. 

The circuit of Fig. 4 is the natural one tq use when 
large temperature variations are expected. An analysis 
will be made of the temperature compensation problem 
for that circuit. For any circuit of this type, the 
voltage across the potentiometer is given by the formula 


c= {[R:/(Ri + ™ [R3/(Rs + Ri) 
If the resistances R, are given increments AR,, then to a 
first approximation the voltage will change by 
de — RAR, RAR; — 
€ (Ri + (Rs + Ry)? 


If the resistances R; are in one piece and have the 
same temperature coefficient, then the change of re- 


sistance due to temperature is given by AR; = aR,;AT 
and the change in meter reading is 


de aR :R,AT—aR,R,AT aR;R,AT — aR;R,AT 
=0 
€ (Ri + R:)? (Rs + 


so that, under these assumptions, temperature com- 
pensation would be perfect. However, in practice it is 
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FIG. 6 
EFFECT OF TEMPERATURE 
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necessary to consider’the effect of the leads which usu- 
ally have a much different temperature coefficient 
from the gages. We denote the lead resistances by R,’ 
and their temperature coefficient by a’ then 


de AT(a — a’)(Re’Ri — Ri'Ry) 
RY Re + 
AT(a — a’)(Ri'Rs — 
(Rs + Rs’ + Ra + Ri’)? 


If we assume that the gages have approximately the 
same resistance R and that R,’ is small compared 
to R, then, neglecting higher order terms in R,’ and 
(Ri — R) 


de Ri’ — Ri’ — Ry + Ri 


a) 


This shows that if the leads are approximately equal, 
the temperature compensation will be satisfactory. 
With only moderate care, the zero drift for 100°C. 
temperature change can be held to an amount corre- 
sponding to Aa/a:-= 1 X 10-5, and by trial-and-error 
adjustment of the leads values of one-fifth of this may 
be reached. Shunting resistance to give a suitable zero 
on the meter has the effect of changing the lead resist- 
ance and should be avoided. The zero can be con- 
trolled by selecting the gages or by adjusting the range 
of the potentiometer. 

Since the sensitivity of a strain gage system depends 
on the applied voltage, it is desirable to know how 
much voltage can be safely used. This is a complicated 
problem since it depends considerably on the heat trans- 


fer characteristics of the installation. Also, the zero 
drift will depend on the temperature compensation and 
the uniformity of the gage temperature. Therefore, 
only experience can supply the final answer to this 
problem. 


Fig. 6 shows the results of some tests. The beam had 
approximately 4 sq.in. effective heat transfer area. 
The area around the gages was wrapped with felt. With- 
out this wrapping readings were extremely unstable, 
even at low voltages. The gages tested were two com- 
mercial Bakelite gages of different sizes and one home- 
made gage cemented with celluloid cement. In Fig. 6, 
e/Le (which is proportional to s/E) is graphed against 
e?/RZD (which is proportional to the energy input per 
unit wire surface area). This choice of abscissa gives 
lines of about the same slope and this suggests that the 
wire temperature is the parameter affecting strain sen- 
sitivity. The breakdown point is entirely different 
on this scalé but seems more nearly to depend on total 
energy input and therefore must be related to beam 
temperature. The breakdown point is marked by pro- 
nounced drifting of the reading under load, probably 
due to slipping. 


In Fig. 7 the results of some tests of the effect of tem- 
perature are given. The value e/Le (which is propor- 
tional to S/E) is graphed against temperature. One of 
the gages was a commercial Bakelite gage, one was a 
commercial celluloid gage, and the third, manufactured 
locally, used celluloid cement. The maximum tempera- 
ture reached with the locally made celluloid gage 
(140°C.) was considerably above that quoted for similar 
commercial gages. 
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three quantities depends on the Mach Number at which the heat 
isadded. Eqs. (a) to (c) were used to obtain Fig.1. Fig. 1 shows 
that at a Mach Number of zero all of the heat added to the fluid 
appears as additional enthalpy (@U + dW). As M is increased, 
alarger and larger proportion of the heat energy goes into kinetic 
energy of directed motion. At M = 1/+/y, all of the heat energy 
goes into kinetic energy of directed motion. For1 > M > 1/ Vy, 
dE/dQ is greater than unity. In other words, all of the added 
heat energy and part of the internal energy of the fluid goes into 
kinetic energy of motion along the pipe. 

(2) It is interesting to note that the condition for the vanish- 
ing of dT, which is the yM? = 1, can be manipulated as follows: 


yM? = y(u?/a*?) = u?/RT = 1 


Therefore, 


u?/2 = RT/2 = '/;(C?/2) 
where C? is the mean square molecular speed of random (thermal) 
motion. Therefore, when y/? = 1, the kinetic energy of the 
directed translatory motion along the pipe of 1 Ib. of the 
gas is equal to the kinetic energy associated (under the assumption 
of equipartition of thermal energy) with each degree of freedom 
in translation of the molecules in 1 Ib. of gas. 

(3) The case of isothermal flow can also appropriately be dis- 
cussed here. Binder has shown? that, for flow with friction and 
heat transfer, the limiting Mach Number for isothermal flow 
is M = 1/ V¥- The physical reason for this limit can be derived. 
The method will only be outlined and the result given briefly. 

It should be remembered that in the discussion in paragraph (1) 
of the present letter the effects of friction were neglected. The 
effects of friction are taken into account by adding a friction term 
to Eq. (4). The specific case of heat transfer by forced convec- 
tion is treated by setting 


mdQ = h(Ty — T,)dS 


where h = heat-transfer coefficient, 7, = temperature of wall, 
T,; = stagnation temperature of gas, and S = surface area of wall. 
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(Continued from page 24) 


EDITOR 


The quantity h can then be eliminated in terms of the friction 
factor f by means of the relation 


= f/2 


It can then be shown that the condition for isothermal flow, 
dT = 0, is that 


— 


T, Di + — 

This equation was employed to obtain Fig. 2, which shows 
T/T; as a function of M for constant temperature. As M ap- 
proaches 1/+/7 from below, the value of 7,,/7; required for iso- 
thermal flow approaches infinity. For M > 1/1/y, negative 
values of 7,,/T; are required. Inasmuch as Ty, and 7; are here 
based on the absolute temperature scale, negative values are 
physically impossible. They are shown by the dashed curve in 
Fig. 2. 

What occurs is this. For all subsonic Mach Numbers, friction 
causes the temperature to decrease. For Mach Numbers up to 
1/+/7, the addition of heat to the fluid causes the temperature to 
increase. It is possible for the two effects to cancel each other 
and for the temperature to remain constant. At M = 1/+/y 
the addition or the subtraction of heat does not affect the te:n- 
perature; therefore, the temperature must decrease because of 
friction. For 1 > M > 1/+/y, the extraction of heat causes the 
temperature to rise, but heat cannot be extracted by forced 
convection at a rate that is great enough to cancel the effect of 
friction. Therefore the temperature cannot remain constant. 

(4) There is a misstatement on page 539, column 2, section 
3(b), of the reference paper which should be corrected. The cor- 
rect statement is: In the supersonic case, when heat is added, 
u decreases and p and p increase. 
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Fic. 2. Ratio of wall temperature to gas stagnation tempera- 
ture as a function of Mach Number for isothermal flow. y = 
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Services 
of the Libraries of the 
Institute of the Aeronautical Sciences 


The services of the Libraries are available to all mem- 
bers of the Institute, to Corporate Members, to advertisers 
in the AzronauticaL ENGingerInG Review and Agro- 
NAUTICAL ENGINEERING CaTALoG, and, under usual library 
limitations, to the public. Four specialized services are 
available. 


The Paul Kollsman Lending Library 


This lending library service makes available, without 
charge, the latest and more important aeronautical books. 

Members may request the loan of any aeronautical or 
technical book they wish to borrow. Through an exchange 
agreement with the Engineering Societies Library, any 
book on general engineering may be borrowed from its great 
collection of over 160,000 volumes. 

A photostating service is available at usual library rates. 

Applications for membership in the library and further 
information will be sent on request. 


The W. A. M. Burden Reference Library 


This reference library contains over 12,000 aeronautical 
books, magazines, pamphlets, and reports gathered from 
world-wide sources and is one of the most complete aero- 
nautical libraries in the world. Material from this library 
fs not available for loan but may be used for reference pur- 


poses. 


The Pacific Aeronautical Library 


6715 Hollywood Boulevard 
Los Angeles 28, California 


Established in cooperation with the aircraft companies 
the library serves. The leading aircraft companies in or 
near Los Angeles participate in its support and operation. 

This service library for aeronautical research is available 
to the public for reading privileges. Source material in- 
cludes aerodynamic and structural research reports, as well 
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2 EAST 64TH STREET, NEW YORK 21, N.Y. 


as books on drafting, production methods, history, and al- 

lied sciences. It furnishes books, periodicals, and pamphlet 
material to the participating aircraft companies to supple 
ment their engineering libraries. 


Technical Information Service 


This service has experienced .personnel under the super- J 
vision of trained aeronautical engineers to compile any in- J 
formation desired. The services range from listing special- 
ized reference books to the preparation of exhaustive 
bibliographies, digesting of reports, and general surveys of J 
any aeronautical subject. Some of the available services § 
are: 


Bibliographies on any aeronautical subject. 
Reports on any aeronautical] subject. 4 
Digests of aeronautical books, papers, periodicals, and refer-@ 
ences. 
Translations. 
Engineering investigations of special aeronautical] subjects. 
Biographies of individuals engaged in aeronautics. 
Photostats of any aeronautical or general engineering 
material. 
Microfilms made on special order. 
Photographs made from the Institute’s photographic collec-@ 
tion. 
Drawings and tracings made. 


In addition to the services mentioned any commission] 
which comes within the scope of the Service will be ag 
cepted. Special arrangements may be made for work re 
quiring several weeks or months. 

Translators are available for accurate transcriptions of all 
foreign language data. Translations are carefully edited] 
by trained engineers. 

Reproductions of any material in the collections off 
the Institute may be ordered at standard photostat rates. 
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